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1.  Background 


Mississippi  State  University  (MSU),  in  collaboration  with  members  of  the  National  Science 
Foundation  Engineering  Research  Center  for  Computational  Field  Simulation  at  Mississippi 
State  University  (MSU)  has  established  a  collaborative  research  program  funded  by  AFOSR, 
McDonnell  Douglas  Corporation,  Boeing  Company  and  Teledyne  Brown  Engineering  (TBE) 
Corporation. 

The  Broad-based  objective  of  this  research  project  is  to  improve  understanding  of  fluid 
phenomena  associated  with  complex  three-dimensional  full  configurations  under  various 
physical  environments  utilizing  hybrid  techniques.  In  particular,  the  research  is  placed  on: 

•  Automatic  hybrid  (Structured,  Unstructured  and  mixed  cells)/generalized  grid 
generation. 

•  Optimal  criteria  for  the  placement  of  structured,  unstructured  and  mixed  cells 

•  Efficient  data  structures  and  data  flow 

•  Euler  and  Navier-Stokes  flow  solvers  for  hybrid/generalized  grids 

•  Training  of  CFD  applications  specialists  and  technology  transfer 

The  key  element  of  this  project  is  the  technology  transfer  providing  one-on-one  interactions 
between  university  personnel  and  scientists  from  various  aerospace  industries.  This  program 
is  a  three-year  collaborative  effort.  This  report  summarizes  the  accomplishments  and 
algorithms  developed  under  this  broad  based  program. 
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2.  Project  Objectives 


The  broad-based  objective  of  this  research  project  is  to  provide  a  forum  through  which  university 
and  industrial  researchers  can  jointly  pursue  research  and  development  pertinent  to  the  simulation 
of  fluid  flowfields  (steady  and  unsteady)  about  complex  configurations  over  a  wide  range  of 
physical  conditions  utilizing  hybrid  techniques.  Upon  completion,  this  research  will  significantly 
reduce  the  “response  time”  for  addressing  complex  applications  and  will  increase  the  quality  of 
simulation  through  the  improvement  of  numerical  algorithms,  surface  rendering  and  grid  generation 
schemes,  and  the  manner  in  which  graphical  techniques  are  used  to  interpret  computed  solutions. 
Specifically,  research  focused  upon  advancing  the  current  state  of  knowledge  regarding  the 
following  items  is  being  conducted: 

•  Grid  Generation  -  Both  structured  and  unstructured  methodologies  will  be  explored  in 
hybrid  development.  In  particular,  the  algorithms  for  refinement,  decomposition, 
redistribution,  and  remapping  algorithms  developed  by  the  principal  investigator  during 
the  1990-93  collaborative  AFOSR  project  will  be  utilized  in  this  development.  Software 
allowing  generation  of  both  structured  and  unstructured  blocks  within  a  given  domain 
decomposition  will  be  developed.  Automization,  optimal  criteria  for 
structured-unstructured  blocking,  and  efficient  data  structures  will  be  explored  in  this 
hybrid  multiblock  development.  The  second  objective  in  the  grid  generation  will  be  to 
generate  hybrid  grids  of  arbitrary  element  types  connected  in  an  unstructured  manner 
(now  are  referred  as  Generalized  Grids).  The  approach  here  will  be  to  utilize  a 
combination  of  the  structured-unstructured  grid  generation  techniques  to  generate 
optimal  grids  for  a  variety  of  complex  configurations. 

•  How  Simulation  Algorithm  -  There  are  numerous  flow  solution  algorithms  for 
structured  and  unstructured  grids  currently  in  widespread  use  around  the  technical 
community.  The  objective  here  will  be  to  develop  a  finite  volume  algorithm  to  solve  the 
Navier-Stokes  equations  on  generalized  grids.  The  approach  will  be  to  utilize  well 
proven  structured  and  unstructured  techniques  and  adapt  them  for  use  with  generalized 
grids.  This  development  will  also  include  an  automatic  grid  adaptation  algorithm  for 
generalized  grids. 

•  Training  of  CFD  Application  Specialists  -  An  often  overlooked  aspect  regarding 
“hardcore”  CFD  applications  is  in  regard  to  what  kind  and  how  much  training  the 
industrial  based  “CFD  Application  Specialist”  needs. 

An  extremely  important  aspect  of  the  present  research  is  in  regard  to  training  of  both  university  and 
industrial  personnel  in  effectively  utilizing  newly  developed  CFD  related  software  using  hybrid 
techniques  for  solving  fluid  flow  problems  associated  with  complex  configurations.  From  a 
day-to-day,  working  level,  people  point  of  view,  it  is  anticipated  that  this  training  will  tend  to  occur 
naturally  as  a  by-product  of  the  collaborative  nature  of  the  efforts  proposed  herein.  However, 
because  it  is  felt  that  training  issues  associated  with  real  world  CFD  applications  are  as  important 
as  the  research  itself,  and  because  this  embraces  of  the  spirit  of  the  Air  Force  initiative  which 
spawned  the  present  proposal,  and  objective  of  this  effort  is  to  take  whatever  steps  deemed  necessary 
to  insure  that  substantial  exchange  of  information  and  personnel  between  this  institution  and  its 
collaborative  partners  take  place. 
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With  the  advent  and  rapid  development  of  supercomputers  and  high  performance  workstations,  it 
is  becoming  possible  to  treat  physical  field  problems  of  engineering  importance  on  very  complex 
regions  by  the  numerical  solution  of  systems  of  partial  differential  equations.  Computational  fluid 
dynamics  (CFD),  for  instance,  has  progressed  to  the  point  where  structured  or  unstructured  flow 
solutions  about  complete  aircraft  configurations  are  possible. 

CFD  involves: 

•  Grid  Generation  -  Creation  of  a  “computational  grid,”  where  the  spatial  region  of  interest 
is  split  up,  or  discretized,  into  many  smaller  regions.  Two  basic  strategies  multiblock 
(composite,  zonal)  (structured  and  unstructured)  are  widely  utilized.  A  multitude  of 
techniques  and  codes  have  been  developed  for  both  structured  and  unstructured  grid 
generation.  However,  grid  generation  is  by  far  the  most  time  (labor)  intensive  part  of  the 
overall  CFD  simulation  process. 

•  Solution  Algorithms  -  Solving  a  set  of  partial  differential  equations  pertinent  to  fluid 
phenomenon  being  analyzed.  Finite  difference,  finite  volume  and  finite  element 
algorithms  with  both  explicit  and  implicit  formulations  have  all  been  used  successfully. 
Codes  utilizing  structured  multiblock  grids  are  widely  used  for  Navier  Stokes  equations 
and  they  have  been  validated  for  a  variety  of  physical  conditions  and  configurations.  In 
the  last  few  years,  unstructured  flow  solvers  have  made  considerable  progress  especially 
in  the  simulation  of  inviscid  (Euler)  flowfields.  The  applicability  of  unstructured 
techniques  to  viscous  and  turbulent  flow  simulations  is  being  researched  by  various 
scientists/engineers. 

•  Numerical  Flow  Visualization  -  post-processing  software  designed  to  allow  the 
rendering  of  a  large  amount  of  digital  information  into  graphical  or  analog  descriptions 
of  a  particular  flowfield  phenomena.  The  widely  used  software  systems  for  post 
processing  are  FAST  (Ref.l)  and  PLOT3D  (Ref.  2).  An  unstructured  capability  has 
recently  been  incorporated  within  FAST.  However,  it  does  not  support  hybrid  grids  and 
further  development  is  needed. 

In  the  spirit  of  AFOSR’s  initiative  to  promote  university-industry  collaborative  research,  a  project 
entitled,  “Numerical  flow  simulations  around  complete  configuration,”  was  executed  during 
1990-93.  The  broad-based  objective  of  this  research  project  was  to  provide  a  forum  through  which 
university  and  industrial  researchers  can  jointly  pursue  research  and  development  pertinent  to  the 
simulation  of  fluid  flowfields  about  complete  flight  vehicles  over  a  wide  range  of  flight  conditions. 
This  collaborative  research  effort  funded  by  AFOSR  and  participating  industries  (McDonnell 
Aircraft  Company,  McDonnell  Douglas  Research  and  Development  and  Teledyne  Brown 
Engineering  Company)  has  been  very  fruitful  for  MSU  and  the  industry  at  Mississippi  State 
University.  Significant  progress  in  the  basic  development  of  surface  grid  technology,  grid  adaptive 
methodology  and  domain  decomposition  techniques  was  realized. 

The  collaborative  project  entitled,  “Hybrid  CFD  techniques,”  was  initiated  as  a  three  year  project 
(1993-96) 
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3.  Technical  Introduction:  Hybrid/Generalized  Grid- 

Flow  Solvers 

Traditionally,  both  experimental  and  theoretical  methods  have  been  used  to  develop 
designs,  and  analysis  for  equipment  and  vehicles  involving  fluid  flow  and  heat  transfer.  With 
the  advent  of  the  digital  computers,  a  third  method,  the  numerical  approach,  has  emerged. 
Although  experimentation  continues  to  be  important,  especially  when  the  flows  involved  are 
very  complex,  the  trend  is  clearly  toward  greater  reliance  on  computer  based  predictions  in 
practical  engineering  applications.  The  development  of  the  high-speed  digital  computers  has 
had  a  great  impact  on  the  way  in  which  principles  from  the  sciences  of  fluid  dynamics  and 
heat  transfer  are  applied  to  problems  of  design,  analysis,  and  manufacture  in  modem  engineer¬ 
ing  practices. 

Many  important  physical  processes  in  nature  are  governed  by  partial  differential  equa¬ 
tions  which  are  usually  difficult  to  solve  or  possess  no  analytical  solutions.  With  the  help  of 
digital  computers,  these  problems  can  be  solved  numerically,  and  the  physical  phenomena  can 
be  simulated  by  computers.  The  numerical  simulations  of  these  equations  are  referred  to  as 
Computational  Field  Simulations  (CFS).  The  CFS  process  involves  numerical  grid  generation 
(generation  of  a  discrete  representation  of  surfaces  or  volumes  for  the  physical  domain),  nu¬ 
merical  solution  algorithms  (numerical  solutions  of  governing  equations  of  fluid  dynamics), 
and  scientific  visualization  (interpretation  of  flow  characteristics  of  physical  domains). 

During  the  past  decade,  computational  simulation  of  fluid  flow  over  complex  configu¬ 
rations  has  progressed  significantly,  and  numerous  notable  successes  have  been  reported  in  the 
literature  (Lohner  [1],  Marcum  and  Weatherill  [2],  Mavriplis  and  Jameson  [3],  Mavriplis  and 
Venkatakrishnan  [4],  Morgon  et  al.  [5],  Venkatakrishnan  and  Mavriplis  [6],  Whitaker  et  al. 
[7],  Whitfield  et  al.[8]) .  However,  the  generation  of  a  high  quality  mesh  for  such  problems 
has  often  been  reported  as  a  pacing  item.  Hence,  much  effort  has  been  expended  to  speedup 
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this  portion  of  the  simulation  process,  resulting  in  several  approaches  to  grid  generation.  The 
generation  of  high  quality  meshes  for  such  complex  configuration  problems  is  the  vital  part  of 
computation  field  simulations  because  the  meshes  are  required  to  accurately  model  the  physi¬ 
cal  domains  and  also  the  mesh  concentration  is  desired  in  the  regions  of  the  domains  where 
flow  features  develop.  The  quality  of  the  mesh  for  a  practical  problem  can  greatly  affect  the 
accuracy  and  efficiency  of  the  solutions  of  the  problem.  The  solution  of  the  governing  equa¬ 
tions  can  be  greatly  simplified  and  the  computational  efficiency  and  accuracy  can  be  improved 
by  well-constructed  set  of  points.  Usually,  the  grid  generation  is  the  most  labor  and  time  con¬ 
suming  part  among  the  whole  computational  field  simulation  processes.  Two  of  the  most 
common  approaches  for  grid  generation  are  based  on  structured  multi-block  (Dannenhoffer 
[9],  Shih  and  Soni  [10],  Thompson  [11])  and  unstructured  procedures  (Kallinderis  [12],  Kha- 
waja  [13],  Lohner  and  Parikh  [14],  Marcum  and  Weatherill  [15],  Mavriplis  [16],  Mavriplis 
[17],  Weatherill  [18]). 

In  the  case  of  structured  grids,  the  physical  domain  of  interest  is  decomposed  into  a 
number  of  quadrilaterals  for  two  dimensions  and  hexahedrons  for  three  dimensions.  These 
regions,  called  cells,  are  numbered  in  such  a  way  that  the  neighbors  of  a  particular  cell  can  be 
determined  in  a  trivial  manner.  The  structured  grid  generation  technology  has  two  major  ap¬ 
proaches:  algebraic  interpolation  approach  which  generates  grids  directly  by  interpolation 
with  the  features  of  economical,  fast,  and  precise  spacing  control,  and  partial  differential 
equations  approach  which  generates  grids  indirectly  by  solving  a  set  of  partial  differential 
equations.  However  the  elliptic  type  Laplacian  partial  differential  equations  approach  can 
avoid  grid  line  crossing  and  the  elliptic  type  Poisson  partial  differential  equations  can  achieve 
grid  line  orthogonality  with  spacing  controls. 

In  the  case  of  unstructured  grids,  the  flow  domain  is  divided  into  triangular  cells  for 
two  dimensional  cases  and  tetrahedrons  for  three  dimensional  cases.  The  cell  numbering  does 
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not  follow  any  particular  pattern.  Therefore,  the  cells  that  neighbor  a  given  cell  are  not  direct¬ 
ly  deducible.  This  necessitates  the  creation  of  an  explicit  connectivity  table  that  contains  this 
information  and  results  in  an  overhead  for  storing  the  same.  However,  unstructured  grids  of¬ 
fer  greater  geometric  flexibilities.  The  generalized  data  structure  of  an  unstructured  grid  is 
very  useful  in  the  refinement  or  de-refinement  of  the  grid  since  the  data  structure  has  to  be 
changed  only  locally.  This  helps  in  adaptation  of  the  grid  to  the  flow  features  by  adding  more 
points  where  the  gradients  of  the  flow  properties  are  significant  ( Baum,  et  al.  [19],  Marcum 
and  Weatherill  [20],  Mavriplis  [17])  and  removing  points  from  the  regions  where  there  are  no 
flow  features  of  interest.  Therefore,  grid  adaptation  is  easier  in  the  case  of  unstructured  grids 
as  compared  to  moving  or  adding  grid  points  to  a  structured  grid. 

A  class  of  problems  of  interest  to  the  Computational  Fluid  Dynamics  (CFD)  communi¬ 
ty  is  the  flow  simulation  over  dynamically  moving  bodies.  Chimera  or  overset  grids  are  the 
primary  choice  for  flow  simulation  of  this  type  of  problems  when  using  structured  grid  meth¬ 
odology.  In  this  approach,  complex  configurations  are  decomposed  into  simpler  geometric 
entities  and  structured  grids  are  generated  around  them.  The  governing  equations  are  solved 
on  these  individual  grids  and  the  solution  variables  are  transferred  appropriately  between  the 
grids  for  the  next  time  step.  The  disadvantage  of  this  approach  is  that  the  interpolation  be¬ 
tween  the  grids  may  not  always  be  conservative.  If  there  are  discontinuities  in  the  flow  vari¬ 
ables  near  the  overlapping  grid  regions,  then  there  may  be  spurious  oscillations  near  the  grid 
interfaces. 

Flow  solvers  based  on  structured  grids  tend  to  be  computationally  more  efficient  than 
those  based  on  unstructured  grids.  High  aspect  ratio  cells  that  are  necessary  for  the  resolution 
of  the  viscous  layers  can  be  easily  generated  for  structured  grids.  The  state  of  the  art  for  struc¬ 
tured  flow  solvers,  (Cooper  and  Sirbough  [21],  Whitfield  et  al.[8],  Thomas,  et  al.[22]),  can 
handle  this  type  of  high  aspect  ratio  cells.  For  complex  configurations,  the  physical  domain 
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has  to  be  decomposed  into  different  sub-domains  (blocks)  and  the  grid  has  to  be  generated 
separately  for  individual  blocks.  In  many  cases,  grid  lines  do  not  exhibit  continuity  across  the 
block  interfaces.  Even  with  this  relaxation  of  the  continuity  of  the  grid  lines  at  the  block  inter¬ 
face,  this  grid  generation  process  is  time  consuming.  Flow  solvers  supporting  non-contiguous 
interfaces  require  specialized  interpolation  procedures  which,  furthermore,  may  not  ensure 
conservation  of  the  flow  variables  at  the  block  interface. 

Flow  solvers  based  on  an  unstructured  mesh  require  more  CPU  time  per  grid  point. 
Furthermore,  the  number  of  grid  points  necessary  to  resolve  the  boundary  layers  using  nearly 
equilateral  triangles  are  enormous,  resulting  in  significantly  higher  CPU  and  memory  require¬ 
ments.  Recently,  success  has  been  reported  in  generating  high  aspect  ratio  unstructured  grids 
for  viscous  simulations  using  an  advancing  normal  point  placement  strategy  (Marcum  [23], 
Pirzadeh  [24]).  In  this  case,  the  triangles  inside  the  boundary  layer  are  very  skewed.  This 
forces  one  to  avoid  the  use  of  a  cell-centered  finite  volume  approach,  as  the  truncation  error  is 
inversely  proportional  to  the  sine  of  the  minimum  angle  of  the  triangles  that  form  the  mesh. 
Consequently,  the  truncation  error  will  be  very  high  in  the  cells  that  are  in  the  boundary  layer 
region.  For  the  node  centered  schemes,  the  diagonal  edges  do  not  contribute  to  the  effective 
flux  balance.  This  results  in  unnecessary  computation  of  the  fluxes  crossing  the  diagonal 
edges. 

Hybrid  or  generalized  element  grid  generation  and  solution  techniques  have  been  de¬ 
veloped  with  the  objective  of  combining  the  attractive  features  of  both  structured  and  unstruc¬ 
tured  techniques.  The  Table  1  compares  the  advantages  and  disadvantages  of  different  grid 
generation  approaches  and  there  is  a  clear  advantage  for  the  hybrid  approach. 
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Structured 

Unstructured 

Chimera 

Hybrid 

Geometric  Flexibility 

— 

+ 

+ 

+ 

Grid  Adaptation 

— 

+ 

— 

High  Aspect  Ratio  Cells 

+ 

— 

+ 

+ 

Moving  Grids 

— 

+ 

+ 

+ 

Interpolation 

— 

+ 

- 

+ 

Memory 

+ 

— 

+ 

? 

CPU 

+ 

— 

+ 

? 

Table  1  Relative  Advantages  of  Different  Grid  Generation  Approaches 


Nakahashi  et  al.  [25]  used  a  zonal  method  for  the  hybrid  grid  generation.  In  their  ap¬ 
proach,  body  fitted  structured  grid  is  kept  in  the  boundary  layer  and  finite  difference  scheme  is 
used  to  solve  the  Reynolds  averaged  Navier-Stokes  equation  in  the  structured  portion  of  the 
grid.  The  rest  of  the  domain  is  filled  with  unstructured  grid  and  the  Euler  equations  are  solved 
in  that  part.  The  interface  between  the  structured  and  unstructured  parts  are  treated  as  explicit 
boundary  condition.  Lohner  [26]  used  a  combination  of  semi-structured  and  unstructured 
grid  for  getting  high  aspect  ratio  cells  in  the  boundary  layer.  The  semi-structured  grid  in  the 
boundary  layer  is  generated  using  the  surface  normals.  In  this  approach  prisms  are  generated 
in  the  boundary  layer  and  are  trimmed  to  avoid  the  grid  crossing.  These  prisms  are  subdivided 
into  tetrahedra  to  only  work  with  tetrahedra  in  the  domain.  The  main  disadvantage  of  the 
method  is  that  the  presence  of  highly  skewed  tetrahedra  in  the  boundary  layer.  Kallinderis  et 
al.  [12],  and  Sharov  and  Nakahashi  [27]  used  a  semi-structured  prisms  in  the  viscous  regions 
and  tetrahedra  in  the  rest  of  the  domain.  In  both  these  cases,  they  used  the  marching  direction 
as  the  surface  normals  to  generate  the  prisms.  Sharov  and  Nakahashi  [27]  used  a  Delaunay 
triangulation  for  the  tetrahedra  generation  while  Kallinderis  et  al.  [12]  used  an  advancing  front 
method. 
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In  the  approach  put  forth  by  Weatherill  [28]  and  Kao  and  Liou  [29],  most  of  the  do¬ 
main  is  filled  with  structured  elements  and  the  different  components  of  the  structured  grid  are 
connected  using  unstructured  grids.  This  involves  the  decomposition  of  the  complex  geome¬ 
tries  into  a  number  of  simple  geometric  entities  followed  by  the  generation  of  structured  grids 
around  them.  One  of  the  structured  grids  is  termed  the  main  grid  and  it  contains  the  remaining 
grids  called  component  grids.  The  main  and  component  grids  are  overlaid  and  a  hole  is  cut  in 
the  main  grid  where  the  structured  -  grid  component  has  to  be  placed.  The  gap  between  the 
main  and  the  component  grids  is  filled  with  an  unstructured  grid.  An  example  of  a  hybrid  grid 
generated  using  this  approach  is  shown  in  Figure  3.1.  As  can  be  seen  from  the  figure,  the 
quality  of  the  grid  at  the  transition  between  the  structured  and  unstructured  grids  is  not  satis¬ 
factory  with  the  area  ratio  of  the  cells  varying  abruptly  across  the  interface.  This  leads  to  a 
higher  truncation  error  during  the  discretization  of  the  governing  equations. 

In  the  present  approach,  structured  grids  are  used  only  near  solid  bodies  and  the  rest  of 
the  domain  is  filled  with  unstructured  grids.  The  structured  grids  are  generated  using  an  ad¬ 
vancing  layer  type  method  and  the  quality  of  the  cells  at  the  interfaces  is  assured  by  checking 
the  aspect  ratios  of  the  cell  (Huang  [30]).  The  hybrid  grids  are  generated  using  a  combination 
of  structured  grid  generator  based  on  an  advancing  layer  method  and  an  unstructured  grid  gen¬ 
erator  based  on  Delaunay  triangulation.  An  example  of  the  hybrid  grid  around  a  two  element 
airfoil  using  the  present  grid  generation  approach  is  given  in  Figure  3.2.  One  can  clearly  ob¬ 
serve  the  improvement  in  grid  quality  in  the  transition  region  between  Figure  3.1  and 
Figure  3.2. 

The  development  of  a  flow  solver  for  generalized  grids  is  a  challenging  problem, 
since  it  has  to  handle  cells  with  arbitrary  number  of  sides.  The  existing  structured  and  un¬ 
structured  grid  algorithms  are  combined  with  the  flow  solver  to  develop  a  flow  simulation 
system  for  hybrid  grids. 
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Figure  3.1  Hybrid  Grid  for  Two  Element  Airfoil 


Figure  3.2  Hybrid  Grid  for  Two  Element  Airfoil 
Using  New  Approach 

In  other  works  related  to  the  flow  solvers  for  hybrid  grids,  the  structured  and  unstruc¬ 
tured  grids  are  treated  as  two  different  blocks  (Mathur  [31],  Tsug  et  al.  [32],  Soetrisno  et  al. 
[33])  and  explicit  boundary  conditions  are  used  to  transfer  information  between  them.  This 


introduces  a  time  lag  between  the  interface  of  the  structured  and  the  unstructured  grids  result¬ 
ing  in  a  possible  convergence  degradation.  In  the  present  work,  the  structured  and  unstruc¬ 
tured  grids  are  treated  as  a  single  block  for  the  purpose  of  flow  computations.  This  is  achieved 
through  a  generalized  data  structure.  Another  existing  hybrid  flow  solver  is  due  to  Parthasara- 
thy  et  al.[34].  In  their  approach,  they  use  Lax-Wendroff  temporal  discretization. 

A  cell  centered  finite  volume  upwind  scheme  based  on  Roe’s  approximate  Riemann 
solver  is  used  for  the  inviscid  flux  evaluation.  The  turbulence  viscosity  is  estimated  using  the 
Spalart-Allmaras  one  equation  turbulence  model  (Spalart  and  Allmaras  [35])  while  Suther¬ 
land’s  law  is  used  for  the  molecular  viscosity.  Both  explicit  and  implicit  schemes  are  imple¬ 
mented  and  validated  with  experimental  data.  The  convergence  of  the  implicit  scheme  using 
approximate  analytical  Jacobians  and  numerical  Jacobians  is  studied.  The  approximate  ana¬ 
lytical  Jacobians  are  calculated  assuming  the  Roe  averaged  matrix  as  constant  and  differentiat¬ 
ing  the  numerical  flux  crossing  the  cell  face.  But  the  numerical  Jacobians  are  evaluated  by 
perturbing  the  the  conserved  variables  by  a  small  amount  and  estimating  the  change  in  the 
numerical  flux  that  cross  the  cell  faces.  For  dynamically  moving  bodies,  the  grid  movement  is 
computed  through  a  spring  analogy  and  the  trajectories  of  these  bodies  are  determined  using 
the  laws  of  classical  mechanics. 
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4.  Hybrid  Grid  Generation 

The  overall  hybrid  grid  generation  procedure  can  be  divided  into  four  basic  steps.  In 
the  first  step  the  complex  geometries  are  decomposed  into  number  of  simples  geometric  enti¬ 
ties.  In  the  next  step  structured  grid  are  generated  around  each  of  these  geometric  entities 
using  an  advancing  layer  type  method.  In  the  third  step,  the  structured  grid  around  all  geomet¬ 
ric  entities  are  kept  together  and  the  overlapping  regions  are  trimmed  by  checking  the  aspect 
ration.  In  the  final  step  the  void  between  the  trimmed  structured  grids  are  filled  using  unstruc¬ 
tured  grids.  These  steps  are  explained  in  detail  in  the  following  sections. 


4.1  Advancing  Front  Structured  Grid  Generator 
An  advancing  front  structured  grid  generator  is  developed  for  generating  a  grid  which 
originally  marches  the  front  out  from  the  boundary  with  appropriate  packing.  This  advancing 
front  structured  grid  generator  can  be  either  used  alone  for  generating  complex  geometry  grid 
efficiently  or  coupled  in  the  hybrid  grid  generation  system. 

Two-dimensional  advancing  front  grid  equations  can  be  written  as  Soni[36]: 


r rj  '  T"t]  —  d2 

or  |r|  x  rn  |  =  A 

or  an  alternative  form  as: 

Xrj  +  y^yrj  =  0 

*2  +  y2  =  d2 

or  x^yn  -  y^xv  =  A 

where 


(4.1) 


(4.2) 
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r  =  (x , y  ) 

(I,  V) 
d: 

A: 


Physical  space 

Computational  space 

Distance  for  marching  out 
Area  for  marching  out 


From  this  system  of  equations,  we  have  already  known  the  object  boundary 
points  distributions  and  the  distance  for  marching  out  in  the  current  level,  that 
is,  and  d  are  known  values.  We  can  solve  the  equations  for  rv  to  locate  the 
points  for  another  level.  By  solving  this  system  of  equations,  the  properties  of  ort¬ 
hogonality  and  point  distributions  in  tj  direction  are  obtained. 

Three-dimensional  front  advancing  grid  equations  can  be  similarly  written  as  Soni 

[36]: 

r(  ■  r(  =  0 

rv  '  r(  ~  0  (4.3) 

=  d 2 

or  as: 

+  zszz  =  0 

Xr,  xK  +  yvy^  +  zvz^  =  0  (4.4) 

x2  +  y^  +  z£  =  d2 

where 

r  =  ( x  ,  y  ,  z  )  Physical  space 
(£,?/,£)  Computational  space 

d  :  Distance  for  marching  out 
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The  equations  described  above  for  two-dimensional  and  three-dimensional  ad¬ 
vancing  front  structured  grid  generator  is  nothing  but  a  set  of  hyperbolic  partial 
differential  equations.  An  alternative  three-dimensional  system  can  be  used  to 
compute  the  high  level  points  more  efficiently. 
r^xrv  =  rt 

(4.5) 

=  dr 

or  (rSx^)  =  v 

y$  Zri  -  y v  *1  = 

Z^XV  ~  ZVX^  =  y ^ 


x^y n  ~  xvy^  =  z% 

x £2  +  y j2  -I-  z £2  =  d 2 

or  xtynzt  +  x%y%zn  +  xvy^z^ 

-  x^ yv -  Xfj y £ z^  -  x^y^zv  =  V 

where 

r  —  ( x  ,  y  ,  z  )  Physical  space 
(  £  ,  tj ,  £  )  Computational  space 

d  :  Distance  for  marching  out 

V  :  Volume  for  marching  out 


(4.6) 


Figure  4.1  and  Figure  4.2  illustrate  how  the  grid  line  or  grid  surface  marches  out  from 
base  line  or  surface  to  another  level. 
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=  constant 


Figure  4.1  Grid  Point  Advances  to  Another  Level  in  Normal  Direction 

for  Two  Dimension  Application. 


,  £  =  constant 

£  =  constant 

Figure  4.2  Grid  Point  Advances  to  Another  Level  in  Normal  Direction 

for  Three  Dimensional  Application. 

4.1.1  Elliptic  Smoothing 

The  advancing  front  structured  grid  generation  has  properties  of  orthogonality  where 
the  grid  points  march  along  the  normal  direction  of  the  previous  level  of  boundary  or  surface, 
and  precise  packing  distribution  which  are  supplied  by  known  distribution  data,  but  does  not 
have  guaranty  of  grid  crossing  resistant.  For  a  smooth  convex  boundary,  there  is  no  problem 
for  generating  a  perfect  grid  by  utilizing  front  advancing  structured  grid  generator.  But  for  a 
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concave  boundary  area,  it  would  cause  grid  line  crossing  after  several  steps  of  grid  line  march¬ 
ing.  To  resolve  this  problem,  an  elliptic  PDE  grid  smoother  is  employed  to  smooth  the  grid  at 
every  new  level  to  prevent  grid  line  crossing.  However,  elliptic  grid  smoother  which  prevents 
grid  line  crossing  diffuses  the  grid  so  that  all  point  concentrations,  and  line  orthogonality  are 
lost  while  solving  the  elliptic  system  iteratively.  Forcing  control  functions  are  introduced  to 
accomplish  field  orthogonality  and  spacing  control  by  elliptic  PDE  grid  system.  Therefore, 
instead  of  solving  the  Laplace  equations,  we  solve  the  Poisson  equations  to  smooth  existed 
grids. 

Equations  of  elliptic  system  in  two-dimensional  and  in  three-dimensional  domains 
can  be  written  as: 


and 


^ocx  "h  £yy  ~  R  (  £  >  V  ) 

Vxx  "b  Vyy  =  Q  (  £  >  V  ) 

( x  ,  y  )  Physical  Space  ( f  ,  rj)  Computational  Space 


(4.7) 


£xx  *b  £ yy  +  £ zz  =  P  (  £  ,  T]  ,  £  ) 

VxX  "b  T]yy  +  T]zz  =  Q  (  £  ,  Tj  ,  £  )  (4.8) 

£xx  -b  £yy  +  Czz  =  P  (  £  >  V  >  £  ) 

( x  ,  y  ,  z  )  Physical  Space  ( £  ,  rj ,  £  )  Computational  Space 


When  P,  Q,  and  R  are  all  zero,  the  above  equations  will  reduce  to  the  Laplace 
equations  which  will  guarantee  the  removal  of  grid  line  crossing,  grid  line  ortho¬ 
gonality  and  a  smooth  grid  in  the  field.  But  the  distributions  and  concentrations 
of  the  grid  will  be  disturbed.  To  preserve  the  grid  distribution,  control  functions 
are  introduced  where  P,  Q,  and  R  are  not  zero. 

The  Poisson  equations  can  be  rewritten  as  a  form  of  P  and  Q  (Soni  [36]), 
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£22  -  2£i2r I,  +  £nrw  =  -  £(r§P  +  r^Q)  (4.9) 

or  a  form  of  <p  and  ip  : 

£22  (r!£  —  —  2£i2r^  +  £11  ( rw  ~  V*  rv )  =  0  (4.10) 


or  an  alternative  form: 

£22  (rg  '  r§  -  ^r|  '  r|>  -  2£l2r!>7  *  r| 

+  £11  (rw  '  r£  “  Vrn  '  rO  =  0  <'4'11') 

£22  ( r£|  '  rtj  ~  <Pr^  '  rt])  ~  2 £12 r|??  ‘  rv 
+  gn(rrm  '  rv  -  ip  rv  •  rv)  =  0  (4.12) 


where 


£11  =  r|  *  r|  » 


£12  =  r£  ‘  rv  and  £22  =  rri  ‘  rv 


r  =  (x,y) 

P  =  —  (p  and 
£  =  I  I  rf  x  r9  |  | 2 


Q  =  - 


£11 

£ 


V» 


Thomas-Middlecuff  proposed  equation  (4.13)  for  the  control  functions  <p  and  ip  by 
assuming  that  g12  equals  to  zero  for  orthogonality  and  one  dimensional  approach  for  each 
direction.  Soni  proposed  equations  (4.14)  for  the  control  functions  <p  and  ip  which  take  con¬ 
tributions  from  arc  length  factor  for  each  direction.  Equations  (4.13)  and  (4.14)  reveal  the  for¬ 
mula  of  the  control  functions  <p  and  ip  respectively  from  Thomas-Middlecuff  and  Soni  [37]. 

<p  =  -=P  and  ip  =  -jp-  (4.13) 

rv 
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(4.14) 


rtt  *  r* 

<p  =  -= — ; — — —  and  ip  = 


r ijr/  ’  Trj 

rv  •  rv 


Soni  proposed  another  form  of  control  functions  which  can  be  derived  from  equations 
(4.15)  and  (4.16)  as  following  (Soni  [38]) 


j,  N  n  (^n )  y 

S22  (  2  0£ll^  2g12  2 

+  £11  (  (&12  \  -  (gf -  V>£l2)  =  0 


(4.15) 


.  ,,,  \  (SllH  ^  _  >.  (^22)^ 

^22^^12^  2  &  §12  '  2 

+  gn((g|2 -  v(g 22^  =  0 


(4.16) 


These  equations  can  take  the  contributions  both  from  the  arc  length  distribution 
factor  and  curvature  factor  into  the  control  functions  with  same  assumptions  for 
orthogonality  and  one  dimensional  approach,  which  are  : 


,  -  rJl  '  rJL  +  Trm  '  r-l 

Y  r%  '  r£  rV  rV 


,  -  rJH  :  rs  JL  rJi  r-l 

v  rv  •  rv  +  ■  ri 


(4.17) 


Similarly,  three  dimensional  0  and  ip  form  of  Poisson  equations  with  assumption  of 
orthogonality  for  grid  lines  inside  can  be  written  as: 


£22  £33  ( r||  0  rl )  +  £11  £33  (ryy  0  rv ) 

+  £11^22  -  0r()  =  0 


(4.18) 


Control  functions  derived  from  equation  (4.18)  are 
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JUji 

JUtJl 

II 

•  r£ 

+ 

rm 

*  r£ 

+ 

rZ 

•  r! 

■ 

rv 

rn 

rc  • 

rll 

rv 

+ 

rm 

rn 

+ 

rv 

^  ' 

'  rZ 

rv 

•  rn 

r£  ' 

-  rZ 

-1- 

rm 

’  r£ 

+ 

' 

' 

r! 

rv  ■ 

rv 

rK 

(4.19) 


One  of  the  characteristics  of  advancing  front  structured  grid  generator  is  that  the  grid 
may  not  look  good  except  for  the  first  grid  line  because  all  other  gird  lines  are  derived  from 
the  previous  line  with  normal  vectors  computed  in  the  previous  line.  So,  when  we  compute  the 
forcing  control  functions  from  the  whole  grid,  it  may  result  bad  control  functions  if  the  grid 
lines  are  crossing  in  the  higher  portion  of  the  grid.  We  can  rewrite  equations  (4.17)  and  (4.19) 
to  prevent  the  problem  by  computing  the  control  functions  from  the  first  grid  line  in  £  direc¬ 
tion  and  from  distribution  in  rj  direction.  Equations  (4.20)  and  (4.21)  are  presented  for  two- 
dimensional  and  three-dimensional  control  functions  respectively.  The  physical  meanings  of 
terms  gn  ,  g2 2  ,  and  g33  are  distances  of  neighboring  points  in  £  ,rj  ,  and  £  directions  re¬ 
spectively.  For  hyperbolic  grid  generation,  the  grid  lines  marching  out  from  the  original  front 
in  rj  direction  in  two  dimensional  generation,  and  in  £  direction  in  three  dimensional  genera¬ 
tion.  We  have  already  known  g2 2  >  and  g33  from  the  user  supplied  marching  distribution  data. 

Thus,  we  can  get  precise  distribution  control  by  solving  Poisson  equations  with  forcing  con¬ 
trol  functions  supplied  in  equations  (4.20)  and  (4.21)  for  two-dimensional  and  three-dimen¬ 
sional  grid  respectively. 
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(4.20) 
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(§nh 

(§33)$ 

§11 
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^33^ 

(§11)?? 

§22 

£33 

§11 

CO 
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1 

(§22^ 

§33 

*n 

§22 

(4.21) 


For  a  free  boundary  condition  type  grid,  the  marching  distances  at  every  grid  points 
are  the  same  according  to  the  global  distribution  data.  That  means  (g22^  °f  equation  (4.20) 
and  (§33)^  ,  and  (g33)v  of  equation  (4.21)  are  zero.  Both  the  equations  could  be  rewritten  as 
equations  (4.22)  and  (4.23). 

1 


0  = 


2  §n 


_  I  (<*28>S 

V  2  §22 


(§u)y_ 

§11 


(4.22) 
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0 

0 
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2  §11 
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2  1  £22 
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(4.23) 
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All  control  functions  are  computed  from  the  original  grid  which  may  have  grid  line 
crossing  somewhere  in  the  grid.  To  ensure  obtaining  good  forcing  control  functions,  0  for  two 
dimensional  grid,  and  0  and  xp  for  three  dimensional  grid  are  computed  at  every  level,  but 
0  for  two  dimensional  grid,  and  6  for  three  dimensional  grid  are  only  computed  at  the  first 
level. 

Smoothing  the  grids  by  solving  the  Poisson  equations  preserves  orthogonality  for  grid 
lines  inside  the  grids,  and  precise  grid  line  concentration,  which  resolves  the  problem  that  con¬ 
centration  of  the  grid  lines  would  be  lost  while  smoothing  the  grids  by  solving  the  Laplacian 
equations. 

Figure  4.3,  Figure  4.4,  and  Figure  4.5  illustrate  the  results  of  front  advancing  grid  gen¬ 
eration  with  elliptic  smoothing  and  without  elliptic  smoothing.  Also,  phenomena  of  elliptic 
smoothing  with  forcing  control  functions  and  without  forcing  control  functions  are  shown 
from  Figure  4.6  to  Figure  4.13. 


(a)  Without  Elliptic  Smoothing.  (b)  With  Elliptic  Smoothing. 

Figure  4.3  Effect  of  Local  Elliptic  Solver  on  Front  Advancing  Grid 

Generation. 
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(a)  Without  Elliptic  Smoothing. 


(b)  With  Elliptic  Smoothing. 


Figure  4.4  Effect  of  Local  Elliptic  Solver  on  Front  Advancing  Grid 

Generation. 


(a)  Without  Elliptic  Smoothing.  (b)  With  Elliptic  Smoothing. 

Figure  4.5  Effect  of  Local  Elliptic  Solver  on  Front  Advancing  Grid 

Generation. 


Figure  4.7  Rocket  Grid  With  Non-Zero  Forcing  Functions 
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4.1.2  Boundary  Conditions 

Cut  line  boundary  points  are  smoothed  by  two  different  approaches  which  are  based  on 
Bezier  curve  and  liner  extrapolation.  Bezier  curve  is  used  for  the  boundary  where  the  first 
point  and  the  last  point  of  the  same  grid  line  overlap  together.  Linear  extrapolation  is  used  to 
extrapolate  the  points  from  inside  the  grid  to  the  boundary.  Without  specifying  the  boundary 
conditions,  the  concave  boundary  of  the  grid  could  cause  the  grid  line  crossing  and  ruin  the 
whole  volume  grid. 

4. 1.2.1  Bezier  curve 

A  cubic  polynomial  is  written  as  parametric  form  in  equation  (4.24).  If  we  have  the 
boundary  conditions  of  a  curve  as  in  equation  (4.25),  we  can  rewrite  the  equation  (4.24)  to 
equation  (4.26)  which  is  Bezier  curve  of  3rd  degree.  By  controlling  the  boundary  conditions 
of  the  curve,  we  can  easily  smooth  the  points  inside  the  curve  for  any  given  location  controlled 
by  parameter  t. 

r(t)  =  a  +  bt  +  ct2  +  dt 3  (4.24) 

r(0)  =  Cx 

r( 1)  =  C4 

(4.25) 

r'(  0)  =  3  (C2  ~  Ct) 
r  ( 1 )  =  3  (  C4  -  C8  ) 

r  ( t )  =  <p1C1  +  02  ^2  +  03  ^3  +  04  ^4  (4.26) 

where 
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4.1.2.2  Linear  Extrapolation 

Linear  extrapolation  is  performed  by  solving  the  intersection  of  a  line  from  inside  the 
grid  and  a  boundary  surface.  The  line  and  surface  can  be  written  as  equations  (4.27)  and  (4.28) 
respectively.  By  solving  these  equations,  we  can  get  the  point,  r,  on  the  boundary  surface. 

r  -  r1  =  mV  (4.27) 

(r  -  r0)  •  e  =  0  (4.28) 

where 

r0  :  known  point  on  the  boundary  surface 

:  known  point  on  the  grid  line 

e  :  normal  vector  on  the  surface 

V  :  vector  of  the  grid  line 

m  :  unknown  parameter  of  line  equation 

4.2  Grid  Trimming 

Structured  grids  around  the  bodies  are  generated  independently.  All  individual  struc¬ 
tured  grids  are  then  supplied  to  the  hybrid  grid  generation  system  for  trimming  the  overlapped 
structured  grid.  Some  overlapped  structured  grids  are  shown  in  Figure  4.16  through 
Figure  4.18.  Figure  4.16  shows  that  six  cylinder  grids  are  generated  and  overlap  with  each 
others.  Figure  4.17  and  Figure  4.18  show  the  three  element  airfoil  structured  grids  and  two 
element  airfoil  structured  grids  respectively. 

Grids  might  overlap  with  others  without  the  control  of  marching  distances  of  each  in¬ 
dividual  process.  So  overlapping  cells  removal  is  the  most  important  concept  of  this  grid  trim¬ 
ming  work.  For  an  optimal  interface  between  quadlateral  cells  and  triangular  cells,  aspect  ra- 
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Figure  4.16  Six  Cylinder  Grids  Overlap  With  Each  Other. 


Figure  4.17  Three  Element  Airfoil  Grids  Overlap  With  Each  Other. 

tios  of  quadlateral  cells  are  also  checked.  The  three  different  steps  in  trimming  the  overlapped 
structured  grids  are  the  following.  First,  overlapping  cells  from  different  sets  of  structured 
grids  are  removed.  Second,  aspect  ratios  of  remaining  cells  are  checked.  Third,  boundary 


Figure  4.18  Two  Element  Airfoil  Grids  Overlap  With  Each  Other, 
points  are  enriched  if  the  distances  between  two  boundaries  are  not  enough  for  optimal  aspect 
ratio  cells  wrapping  around  the  grids. 

Overlapping  cells  removal  process  is  based  on  the  comparison  of  distances  between 
nodes  of  cells  from  different  set  of  structured  grids  and  solid  boundaries.  This  process  is  the 
initial  step  of  grid  trimming  work,  which  only  ensures  that  all  overlapping  cells  are  removed 
and  all  individual  grids  are  separated  from  others  by  a  clear  area,  or  a  void  area.  An  aspect 
ratio  checking  process  then  performs  to  remove  all  cells  which  does  have  an  aspect  ratio  less 
than  unity.  High  aspect  ratio  cells  are  kept  around  individual  solid  boundaries  for  efficient 
viscous  and  turbulent  flow  simulations  (Marcum  [39]).  Low  aspect  ratio  cells  far  away  from 
the  solid  boundaries  are  removed.  Because  low  aspect  ratio  cells  cause  bad  grid  interface  be¬ 
tween  quadlateral  and  triangular  cells,  which  affect  the  smooth  area  transition  of  the  cell  in 
structured  to  unstructured  region. 

Boundary  points  are  enriched  locally  where  the  distance  between  two  solid  boundaries 
is  not  far  enough  to  have  perfect  aspect  ratio  cells  wrapping  around  the  boundaries  (Marcum 
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[39]).  If  a  region  between  two  boundaries  is  very  narrow,  the  local  flow  solution  is  very  criti¬ 
cal  to  the  global  one.  Figure  4.19  to  Figure  4.22  show  the  comparison  of  the  grid  trimming 
with  and  without  aspect  ratio  checking.  Figure  4.23  and  Figure  4.24  illustrate  the  trimmed 
grids  of  two-element  and  four-element  airfoils.  Boundary  point  enrichment  effect  are  shown 
in  Figure  4.25  to  Figure  4.32. 

The  void  between  the  trimmed  structured  grids  are  filled  with  triangles.  Delaunay 
triangulation  is  used  to  generate  triangular  cells. 


Figure  4.19  Grid  Trimming  Without  Aspect  Ratio  Checking. 

Results  of  two  dimensional  and  three  dimensional  grids  utilizing  front  advancing  and 
elliptic  PDE  smoothing  grid  generator  are  shown  below  from  Figure  4.33  to  Figure  4.40. 

Hybrid  grid  arounf  multiple  cylinders  and  around  twoarbitrary  geometries  hybrid  grids 
are  shown  in  Figure  4.41  and  Figure  4.42  respectively.  Hybrid  grids  around  multi-element 
airfoils  of  practical  importance  to  aerospace  engineering  applications  are  shown  in 
Figure  4.43  to  Figure  4.45. 
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Figure  4.20  Grid  Trimming  With  Aspect  Ratio  Checking. 


Figure  4.21  Grid  Quality  Is  Not  Good  Without  Aspect  Ratio  Checking. 


Figure  4.22  Grid  Quality  Is  Improved  With  Aspect  Ratio  Checking. 


Figure  4.23  Trimmed  Grids  For  Two— Element  Airfoil. 


Figure  4.24  Trimmed  Grids  For  Four— Element  Airfoil. 


Figure  4.25  Grid  Trimming  Without  Boundary  Point  Enrichment. 


Figure  4.26  Grid  Trimming  With  Boundary  Point  Enrichment. 


Figure  4.27  Grid  Trimming  Without  Boundary  Point  Enrichment. 


Figure  4.28  Grid  Trimming  With  Boundary  Point  Enrichment. 


Figure  4.29  Grid  Trimming  Without  Boundary  Point  Enrichment. 
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Figure  4.30  Grid  Trimming  With  Boundary  Point  Enrichment. 


Figure  4.31  Grid  Trimming  Without  Boundary  Point  Enrichment. 


Figure  4.35 


'wo— Element  Airfoil  Grid. 
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Figure  4.37  Three  Dimensional  Surface  Grid 


Figure  4.39  Three  Dimensional  Arbitrary  Geometry  Grid. 
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Figure  4.40  Three  Dimensional  Arbitrary  Geometry  Grid. 


Figure  4.41  Nine  Cylinder  Hybrid  Grid. 


\/J£ 


Figure  4.42  Two  Arbitrary  Geometries  Hybrid  Grid. 


Figure  4.43  Two  Element  Airfoil  Hybrid  Grid. 


Figure  4.44  Three  Element  Airfoil  Hybrid  Grid. 


Figure  4.45  Four  Element  Airfoil  Hybrid  Grid. 


5.  Hybrid  Flow  Solver 

A  hybrid  grid  is  a  combination  of  structured  and  unstructured  grids.  The  number  of 
edges  (2D)  or  faces  (3D)  that  form  a  cell  can  vary  from  one  cell  to  another.  The  arbitrariness 
in  the  number  of  faces  for  a  given  cell  makes  the  finite  volume  method  the  natural  choice  for 
solving  the  governing  equations.  The  finite  volume  schemes  are  based  on  a  discretization  of 
the  integral  form  of  the  equations.  The  governing  equations  used  in  the  present  study  are  the 
Navier-Stokes  equations.  The  following  section  presents  the  Navier-Stokes  equations  in  the 
integral  form. 

5.1  Governing  Equations 

The  non-dimensionalizations  are  done  based  on  freestream  conditions.  The  character¬ 
istic  velocity  is  taken  as  the  freestream  velocity.  Neglecting  the  body  forces,  the  non-dimen¬ 
sional  form  of  the  Navier-Stokes  equations  can  be  written  in  vector  notation  as, 

ft  I  QdQ  +  6  F(Q)  • 

JQ  JdQ 

where  F(Q)  •  nds  =  T  =  [f(Q)i  +  g(Q)j_  +  KQ)k)  '  nds  (5.2) 

F\Q )  •  nds  =  Tv  =  (f\Q)i  +  g\Q)j_  +  h\Q)  k )  •  nds  (5.3) 

n  =  nxi  +  nyj_  +  nzk  (5.4) 

The  vector  n  is  the  outward  unit  normal  to  the  control  surface  with  nx,  ny,  and  nz 
as  the  components  in  x,  y,  and  z  directions  and  ds  is  the  cell  face  area  for  3— dimen¬ 
sions  and  the  edge  length  in  2-dimensions.  The  conserved  variables  Q,  the  con¬ 
vective  flux  vectors/,  g,  and  h  and  the  viscous  flux  vectors  f,  gv,  and  hv  are  de¬ 
fined  as, 


J  dQ 


nds  =  q>  F\Q)  •  nds 


(5.1) 
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where  E  is  the  total  energy  per  unit  volume  and  is  defined  as, 


E  =  pet  =  — j  +  p 


U2  +  V2  +  w2 


(5.5) 


The  grid  speeds  in  the  x,  y,  and  z  co-ordinate  directions  are  oq,  yt,  and  z^  respec¬ 
tively  and  the  variables  U,  V  and  W  are  defined  as,  u  —  xt,  v  —  yt,  and  w  —  zt 
respectively.  The  ratio  of  specific  heats  y  for  standard  air  is  taken  as  1.4. 

The  components  of  viscous  shear  stress  tensor  T,y  represents  the  shear  stress  acting  in  a 
plane  of  constant  i  in  the  direction  of  j,  with  i  and  j  as  either  x,  y,  or  z.  The  heat  transfer  by 
conduction  in  the  i  direction  is  represented  as  q*.  The  components  of  viscous  shear  stress  ten¬ 
sor  r,  and  the  heat  transfer  q  are  defined  as, 
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where  is  the  freestream  Mach  number.  The  turbulence  viscosity  fit  is  calcu¬ 
lated  using  the  Spalart-Allmaras  one  equation  turbulence  model  (  Spalart  and 
Allmaras  [35]  ).  The  turbulent  Prandtl  number  Prt  is  taken  as  0.92. 

The  nondimensional  form  of  the  equation  of  state  for  an  ideal  gas  is  written  as, 

d  =  -£i-  (5.6) 


The  nondimensional  laminar  viscosity  coefficient  is  given  by  the  relation 

(1  +  c3)  P/2 
**  ~  TTT-i 


(5.7) 


where  c3  =  -r~  and  the  reference  temperature  Too  is  taken  as  300  K. 
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5.2  Cell  Face-Based  Data  Structure 


Hybrid  grids  consist  of  polygons  with  an  arbitrary  number  of  sides.  This  makes  it  nec¬ 
essary  to  use  a  generalized  data  structure  to  store  the  grid  information.  The  grid  information 
can  be  stored  in  two  different  ways:  (i)  Cell-based  data  structure  (ii)  Cell  face-based  data 
structure.  These  two  approaches  are  best  explained  by  an  example. 

Consider  a  two-dimensional  grid  that  consists  of  two  cells,  one  with  five  sides  and  the 
other  with  three  sides  (Figure  5.1).  The  nodes  are  denoted  by  Nj,  N2,  ..,iV<5,  the  cells  by  Cj,  C2 
and  the  edges  by  ej,  e2, ej.  In  the  cell-based  data  structure  the  grid  is  represented  as  shown 
in  Table  5.1.  The  nodes  that  form  a  cell  are  usually  ordered  in  either  a  clockwise  or  a  counter¬ 
clockwise  direction.  In  this  case,  the  counter-clockwise  ordering  of  nodes  was  chosen. 


Figure  5.1  Example  of  Hybrid  Grid 
Table  5.1  Cell  Based  Data  Structure 


Cell  No. 

No.  of  Nodes 

Node  1 

Node  2 

Node  3 

Node  4 

Node  5 

Ci 

5 

Nj 

n2 

n3 

n4 

N5 

c2 

3 

n2 

n6 

n3 

In  2-dimensional  cases  the  cell  face-based  data  structure  will  change  to  an  edge-based 
one.  For  the  edge-based  data  structure,  every  edge  has  four  pieces  of  information  associated 
with  it.  They  are  the  nodes  that  form  the  edge  and  the  cell  numbers  on  either  side  of  the  edge. 
For  the  boundary  faces,  the  boundary  condition  information  is  stored  in  place  of  the  cell  num¬ 
ber  on  the  right  side.  The  left  and  right  sides  of  an  edge  (2D)  or  a  face  (3D)  are  defined  in  the 
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following  manner.  When  one  traverses  an  edge  from  the  first  node  to  the  second  node,  the  cell 
on  the  left  hand  side  is  denoted  as  the  left  cell  and  the  cell  on  the  other  side  as  the  right  cell.  In 
three  dimensions,  the  nodes  that  form  the  cell  faces  are  ordered  in  such  a  manner  that  the  num¬ 
bering  is  counter-clockwise  with  respect  to  one  cell  and  clockwise  with  respect  to  the  other. 
The  cell  with  respect  to  which  the  node  numbering  is  counter-clockwise  is  taken  as  the  left 
cell.  In  the  edge-based  data  structure,  the  grid  given  in  Figure  5.1  is  represented  as  shown  in 
Table  5.2. 

Table  5.2  Edge-Based  Data  Structure 


Edge  Number 

First  Node 

Second  Node 

Cell  on  left 

Cell  on  right 

ei 

Nj 

n2 

Ci 

be 

e 2 

n2 

n3 

Cl 

C2 

23 

n3 

n4 

Cl 

be 

e4 

n4 

N5 

Cl 

be 

25 

N5 

Ni 

Cl 

be 

26 

N6 

n3 

C2 

be 

e 7 

n2 

n6 

C2 

be 

In  Table  5.2 ,  be  represents  the  index  for  boundary  conditions. 

The  advantage  of  the  edge  based  system  is  that  it  is  independent  of  the  number  of 
edges/faces  of  the  cells.  Also,  the  integration  over  the  control  surfaces  is  made  easier  by  em¬ 
ploying  this  type  of  data  structure.  This  can  be  better  explained  by  a  method  for  evaluating  the 
area  of  the  cells.  The  area  of  a  cell  can  be  estimated  by  the  simple  integral  as, 

Area  =  <p  xdy  =  (5.8) 

~  i =  1 

where  k  is  the  number  of  cell  edges.  Using  the  edge-based  data  structure,  the 
above  integral  can  be  evaluated  for  the  entire  grid  using  the  following  procedure. 
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Initialize  all  areas  to  zero 
Loop  over  all  edges 

Calculate  the  average  value  ofx  coordinate,  xe 
Calculate  the  change  in  y  coordinate,  dy 
Area  of  cell  on  left  =  Area  of  cell  on  left  +  xe  *  dy 
Area  of  cell  on  right  =  Area  of  cell  on  right  -xe  *  dy 
end  loop. 

The  same  principle  can  be  applied  for  flux  evaluation  too.  Once  the  the  flux  is 
evaluated  for  a  face,  appropriate  contributions  can  be  added  to  the  cells  on  either 
side  of  the  face. 


5.3  Spatial  Discretization 

Finite  volume  ( Barth  [42],  Whitakar  [43] )  and  finite  element  based  schemes  ( Jame¬ 
son  et  al.  [44],  Marcum  and  Weatherill  [2],  Morgon  et  al.  [5]  )  are  the  most  common  ap¬ 
proaches  used  for  solving  the  governing  equations  on  unstructured  grids.  Finite  volume 
schemes  are  best  suited  for  hybrid  grids  because  a  typical  hybrid  grid  is  an  agglomeration  of 
polygons  with  different  number  of  sides  per  polygon.  There  are  basically  two  different  ap¬ 
proaches  for  storing  the  conserved  variables,  one  is  a  node  based  approach  while  the  other  is  a 
cell  based  approach.  Figure  5.2  shows  the  areas  (shaded  region)  that  are  used  for  estimating 
the  averaged  values  of  the  conserved  variables  for  the  two  different  approaches.  In  the  node 
based  method  the  area  considered  for  this  averaging  is  the  area  around  the  nodes.  This  area  is 
termed  dual  area  since  this  is  considered  a  dual  to  the  area  of  the  cell.  For  two  dimensional 
cases  this  area  is  taken  as  the  area  enclosed  by  the  polygon  formed  by  connecting  the  cell  cen¬ 
ters  of  the  neighboring  cells  to  the  midpoint  of  the  corresponding  edges.  In  the  second  ap¬ 
proach,  the  values  stored  at  the  cell  center  are  taken  as  the  cell  averaged  values.  The  area  for 
the  averaging  is  taken  as  that  of  the  cell  itself. 
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In  the  present  study,  the  cell  centered  approach  is  used  for  storing  the  conserved  vari¬ 
ables  and  for  the  subsequent  solution  procedure.  The  cell  averaged  flow  variables  at  the  cell 
center  are  calculated  as  an  integral  average,  i.e., 

Q(Ct)  =  £  Q(x,y,z)  dQ  (5.9) 

1  JQ 

where  Q  is  the  cell  number  and  Vi  is  the  cell  volume. 


(a)  Node  Based  Data  (b)  Cell  Based  Data 

Figure  5.2  Areas  Considered  for  Averaging  the  Conserved  Variable 

in  Node  and  Cell  Based  Schemes 


5.3.1  Discretization  of  Convective  Terms 
The  governing  equations,  equation  (5.1),  without  the  viscous  terms  are, 

QidQi  +  <1  F(Q).nds  =  0  (5.10) 

JQi  J  di 2 

Using  the  cell  centered  finite  volume  approach,  a  semi-discretized  form  of  Equa¬ 
tion  (5.10) is 


(5.11) 
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The  indices  i  andy  denote  the  cell  and  face  numbers  respectively  and  k  is  the  total 
number  of  faces  of  the  ith  cell.  The  right-hand-side  of  equation  (5.11)  represents 
the  flux  balance  for  the  ith  cell. 

The  differential  form  of  the  equation  (5.10)  is  hyperbolic  in  time:  therefore  the  flux 
evaluation  is  based  on  the  direction  of  propagation  of  the  information.  This  is  done  using 
upwinding  based  on  the  eigensystem  of  the  governing  equations. 

The  numerical  flux  crossing  a  cell  face  is  calculated  as  the  exact  solution  of  the 
approximate  Riemann  problem  (  Roe  [45)  ).  The  Riemann  problem  is  defined  as  the  con¬ 
servation  law  together  with  particular  initial  data  consisting  of  two  constant  states  separated 
by  a  single  discontinuity  ( LeVeque  [46] ).  The  basic  principle  behind  the  Roe’s  approximate 
Riemann  solver  is  explained  below  using  the  one  dimensional  hyperbolic  system  of  conserva¬ 
tion  laws, 

Tt  + I  =  0  (6-12) 

with  the  initial  conditions,  as  given  in  Figure  5.3. 


q(x,  0)  = 


qL  x  <  0 
qR  x  >  0 


(5.13) 


4r 


x-0 


Figure  5.3  Initial  Conditions  for  the  Conservation  Laws 


The  above  system  of  conservation  laws  can  be  written  as 


53 


(5.14) 


dq 

dt 


0 


where  a  = 

dq 

Roe’s  approximate  Riemann  solver  finds  the  exact  solution  of  the  following  linear  hy¬ 
perbolic  system. 

f  +  4  =  0  (5.15) 

where  a  is  a  locally  constant  matrix  and  has  to  satisfy  the  following  properties 
(Roe  [45]). 

1)  It  represents  a  linear  mapping  from  the  vector  space  q  to  the  vector  space/ 

2)  AsqL->qR->q,  a(qL>  qR)  ->  a(q) 

3)  For  any  qL,  qR,  f(qR)  -f(qO  =  a(qR,qL)  .(qR-qL) 

4)  The  eigenvectors  of  a  are  linearly  independent. 

The  flux  vector  at  x-0  is  then  given  by 

fx= o  =  |(fL  +  /*)  “  («*  ~  4l)  (5.16) 

In  the  case  of  three  dimensional  Euler  equations  the  Jacobian  matrix  A  is  defined  as 

A  =  and  the  Roe  averaged  variables  appearing  in  A  are  evaluated  using  the  following 
°Q 

expressions  ( Roe  [45] ). 

u 

_  _  iPhvL  +  vR'lpR 

Jp~L+  vS  W 

h  =  v tpjdlk  +  !ieJre. 

JFl  +  K 

where  h  denotes  the  total  enthalpy  per  unit  volume. 
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Using  the  approximate  Riemann  solver,  the  flux  through  a  cell  face  is  given  by 

Zij  =  \[mR)  +  $(QD  ~  \A\(Qr-  Ql)]  (5.17) 

where  I A  I  =  T \A  I  T~l.  T  is  a  matrix  whose  columns  are  the  right  eigenvectors  of 
A,  T~]  is  a  matrix  with  its  rows  as  the  left  eigenvectors  of  A,  and  1/1 1  is  a  diagonal 
matrix  whose  elements  are  the  absolute  values  of  the  eigenvalues  of  5.  The  ma¬ 
trix  A  is  evaluated  using  the  Roe  averaged  variables  defined  earlier. 

The  last  term  of  equation  (5.17)  is  separated  into  three  vectors  based  on  the  three  dis¬ 
tinct  eigenvalues  of  A.  These  vectors  are  given  by 


\A\{Qr-Qj)  =  \AFX  \  +  \AF23A  I  +  \AF5 1 


(5.18) 


\AFl  I 
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cnz 
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+  P 
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\  2  c2 


1_  ' 
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v  +  cny 
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v  j 


where  U  is  the  velocity  normal  to  the  face  which  is  given  by  U  =  unx  +  vny+wnz 
+  nt .  The  symbol  A  represents  the  jump  in  the  values  of  the  flow  variables  across 
the  face  and  is  defined  as  A  ()  =  (  )r-(  )l-  The  quantities  c,  and  q  are  the  speed  of 
sound  and  total  velocity  respectively,  evaluated  using  the  Roe  averaged  variables. 
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5.3.2  Higher  Order  Schemes 


For  a  first  order  accurate  scheme,  the  flow  variables  are  assumed  to  be  constant  within 
each  cell.  The  states  of  the  flow  variables  on  the  left  and  the  right  sides  of  a  cell  face,  for  the 
use  in  the  approximate  Riemann  solver,  are  taken  as  the  values  of  the  cell  averaged  data  on  the 
respective  sides.  This  procedure  avoids  the  creation  of  local  maxima  or  minima  of  the  flow 
variables  and  preserves  monotonicity  ( Barth  [42] ). 

For  second  order  schemes,  the  flow  variables  are  assumed  to  be  distributed  linearly 
within  each  cell  and  the  linear  distribution  is  reconstructed  from  the  cell  averaged  values. 
Figure  5.4  shows  the  distribution  of  the  flow  variables  for  first  and  second  order  schemes  for 
a  one  dimensional  case.  From  the  figure  it  is  evident  that  the  function  can  be  better  repre¬ 
sented  using  linear  reconstruction  than  the  piecewise  constant  representation.  This  in  turn  im¬ 
proves  the  solution  accuracy.  During  the  solution  process  the  flow  variables  at  a  cell  face  are 
extrapolated  from  the  cell  averaged  values  using  a  linear  reconstruction  procedure  (Barth 
[42]). 


(a)  Piecewise  Constant  Data  (b)  Linear  Reconstruction 

Figure  5.4  Distribution  of  Flow  Variables  for  First  and 

Second  order  Scheme. 
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The  piecewise  linear  reconstruction  of  the  conserved  variables  from  the  cell  averaged 
data  is  done  using  the  Taylor’s  series  expansion.  The  Taylor’s  series  expansion  for  a  function 
of  multiple  variables  can  be  written  as, 

Q(x,y,z)  =  Q{xt  ,yitZi )  +  V<2(*;  tyitzj)  -Ar  +  0(  (Ar)2 )  (5.19) 

where  VQ  is  the  gradient  of  Q  and  zlris  the  vector  from  the  center  of  the  cell  ( X[,  y,% 
Zi )  to  the  desired  point  ( x,  y,  z )  and  is  defined  as 

Ar  =  (x  -  xt)  i  +  iy-y^i+iz-  zt)  k  (5.20) 

Using  the  definition  of  Ar  given  in  equation  (5.20),  equation  (5.19)  [neglecting  higher 
order  terms]  can  be  written  as 

Q(x,y,z)  =  Q(x,y,z)  +  ^  (x  -  xj  +  (y  ~  y,)  +  (z  -  zt)  (5.21) 

The  gradient  of  Q  at  the  cell  center  is  estimated  using  Green’s  theorem  (Warsi  [40]). 
Green’s  theorem  relates  the  gradients  within  a  control  volume  to  the  surface  integral  as: 

VQ  =  i  J>  Qnds  (5-22> 

1  JdQ 

The  control  volume  for  integrating  equation  (5.22)  is  taken  as  the  volume  of  the  cell 
itself.  This  necessitates  the  estimation  of  the  properties  at  the  nodes.  The  properties  at  the 
nodes  are  calculated  using  the  weighted  average  of  the  properties  within  the  cells  surrounding 
that  node.  Figure  5.5  shows  the  cells  that  are  considered  for  the  weighted  average  of  the  prop¬ 
erties  at  the  node  Nl.  The  weight  is  taken  as  the  inverse  of  the  distance  of  the  node  from  the 
cell  center  ( Frink  et  al.  [48] ).  The  resulting  equation  for  the  properties  at  the  nodes  are  writ¬ 
ten  as. 
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(5.23) 


Qn\ 


where  Nl  is  the  node  number,  cj  the  neighboring  cell  number,  nj  the  total  number 
of  cells  surrounding  the  node  Nl,  and  rj  the  distance  of  the  node  from  the  cell  cen¬ 
ter,  i.e., 

rj  =  /( *c  -  xj  )2  +  ( -  Jj  )2  +  ( Zc  -  Zj  )2  (5.24) 

where  ( yc,  Zc  )  is  the  cell  centroid  and  (  Xj,  yj,  Zj  )  are  the  co-ordinates  of  the 
node. 


Figure  5.5  Cells  Contributing  to  Node  Nl  for  the 

Weighted  Average 


The  discretized  form  of  the  equation  (5.22),  which  is  used  for  the  gradient  calculation, 
can  be  written  in  terms  of  the  properties  at  the  nodes  as, 

nk 

(V0«  =  y  X  QfR*  dsk  (5.25) 

1  k=  i 
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where  is  the  number  of  cell  faces  for  the  cell  i  and  <2/  is  the  average  value  of  the 
flow  variables  at  the  nodes  which  forms  the  cell  face  k. 

The  derivative  of  the  conserved  variables  in  the  three  co-ordinate  directions  can  be 
explicitly  written  as  (from  equation  (5.25)), 


dQi  = 

dx 

nk 

y  X  Qf(n*>k 

(5.26) 

1  jt=i 

dQi  = 
dy 

y:  X 

(5.27) 

1  *=  i 

dQi  = 
dz 

nk 

y  X  2/0*^* 

1  k= i 

(5.28) 

Using  equations  (5.26)— (5.28)  in  equation  (5.21),  the  linear  distribution  of  the 
properties  with  a  cell  can  be  calculated. 

One  of  the  disadvantage  of  this  approach  is  that  the  information  from  the  wrong  side  of 
the  cell  face  also  contribute  to  the  weighted  averaged  values  of  the  conserved  variables  at  the 
nodes.  This  can  be  overcome  by  using  upwind  biased  gradients  as  discussed  by  Cabello  [49]. 
The  higher  order  calculations  based  on  the  upwind  biased  gradients  are  not  implemented  in  the 
present  study  because  of  the  higher  memory  overheads  required  to  store  the  information  about 
the  cells  used  for  the  upwind  biased  gradient  calculation  and  the  coding  complexity. 

5.3.3  Limiter  Functions 

For  cell  centered  schemes,  the  cell  averaged  values  are  stored  at  the  cell  center.  The 
gradient  within  a  cell  is  calculated  based  on  the  values  in  the  neighboring  cells.  Using  the 
gradients  and  the  cell  averaged  values,  the  linear  distribution  is  reconstructed  as  shown  in 
Figure  5.6.  From  Figure  5.6(a),  one  can  see  that  local  extrema  are  created  at  locations  A  and 
B.  At  location  A  the  reconstructed  value  is  more  than  the  cell  averaged  values  of  the  neighbor¬ 
ing  cells.  Similarly,  at  point  B  the  reconstructed  value  of  the  variable  is  lower  than  that  of  the 
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neighboring  cells.  These  non-physical  values  will  cause  spurious  oscillations  of  the  flow 
variables  near  regions  of  high  solution  gradient.  This  needs  to  be  avoided  to  preserve  mono¬ 
tonicity  and  to  avoid  the  development  of  spurious  oscillations  near  discontinuities.  This  is 
achieved  by  the  application  of  a  limiter  function  (  Barth  [42] ).  However  it  should  be  noted 
that  the  limiter  function  reduces  the  order  of  accuracy  in  those  regions. 


(a)  Before  Applying  Limiter  Function 


Figure  5.6  Effect  of  Limiter  Function  on  Linear  Reconstruction 


The  limiter  function  is  applied  to  the  Taylor’s  series  expansion  to  avoid  the  creation  of 
local  minima  and  maxima.  The  Taylor’s  series  expansion,  given  in  equation  (5.19),  is  modi¬ 
fied  with  the  limiter  function  as, 

Q(x, y, z )  =  Q(x(, yt, zt)  +  <pt  V £>(*/, y,-, Zj)  .Ar+0(  (Ar)2 )  (5.29) 

The  limiting  function  <pi  is  constructed  to  satisfy  the  monotonicity  principle  and 
has  to  satisfy  the  following  condition. 

If  Qimin  =  min(Qci,Qadj) 
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Qmax  -  max  (  Qa  Q^.  ) 

then  Qimin  s=  Q(x,y)  Q™* 

Ci  represents  the  cell  for  which  the  limiter  function  is  being  evaluated  and  the 
subscript  adj  represents  the  cells  adjoining  to  Q.  Two  different  limiter  functions 
were  considered  and  their  effect  on  the  convergence  of  the  numerical  scheme  were 
studied.  The  two  limiter  are  discussed  below. 


5. 3.3.1  Barth’s  Limiter 

A  limiter  function  was  developed  by  Barth  [42]  to  satisfy  the  above  mentioned  proper¬ 
ties  and  is  defined  as 


min^l, 

min^l, 

1 


Qf™  ~  Qy 

Qi  —  Qy 

Qfn  ~  Qy 
Qi  ~  Qij 


)  if  e,  -  Qij  >  o 

j  if  Qi  -  Qn  <  o 
if  Qi  ~  Qv  =  0 


and  =  min(0y) 

In  the  above  relations,  Qij  is  calculated  using  equation  (5.19). 


(5.30) 


(5.31) 


5.3.3.2  Venkatakrishnan’s  Limiter 

Even  though  Barth’s  limiter  gives  accurate  results  for  the  flow  solutions,  the  conver¬ 
gence  stalls  after  a  few  orders  of  magnitude  drop  in  the  residual.  To  avoid  this,  Venkatakrish- 
nan  [50]  proposed  another  limiter  for  unstructured  grids.  This  limiter  is  basically  an  extension 
of  the  van  Albada  limiter,  developed  for  structured  grids.  The  new  limiter  is  defined  as 


<Pij= 


1 


A- 


(A\  +  e2)A-  +  2A2-A  + 
A\  +  2  A2-  +  A- A  +  +  e: 


(5.32) 
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where  s2  is  taken  as  (KAx)3  and  Ax  represents  the  average  edge  size  of  triangles, 
i.e.,  if  the  physical  domain  is  covered  with  same  number  of  equilateral  triangles 
of  equal  size,  then  Ax  is  the  edge  size  of  those  triangles  and  K  is  a  threshold  pa¬ 
rameter.  When  £=0.0,  the  limiter  will  be  active  everywhere  in  the  field,  whereas 
a  very  high  value  of  K  implies  effectively  no  limiting.  Typically  K  is  set  to  0.1  to 
0.3. 

The  other  variables  appearing  in  equation  (5.32)  are  defined  as  follows, 

A  -  =Qij-  Qi  ) 

!  if  Qij-Qi>  0  (5.33) 

A  +  =  Qf™  -  Qi  ) 

A  -  -  Qij  -  Qi  ) 

\  if  Qij  ~  Qi  <  0  (5.34) 

A+=  Q?'n  -Qi  ) 

In  order  to  avoid  division  by  a  small  value  ofzL  in  equation  NO  TAG  it  is  replaced  by 
SIGN(A_)  (\A_  \  +  (o)  and  co  is  set  to  10-12,  where  SIGN(x)  is  the  sign  of*.  Results  are 
presented  for  the  effect  of  these  limiters  on  convergence  in  Section  6.2. 

One  of  the  drawbacks  of  the  Roe  scheme  is  that  it  allows  the  creation  of  nonphysical 
expansion  shocks  in  some  cases.  Another  drawback  ( more  likely  at  higher  Mach  numbers)  is 
that  it  admits  spurious  solutions  along  the  stagnation  line  ahead  of  bow  shocks,  which  is 
known  as  the  carbuncle  phenomenon  (  Quirck  [47] ).  In  order  to  avoid  these  phenomena,  the 
wave  speeds  are  modified  when  the  absolute  value  of  the  eigenvalues  of  A  are  less  than  0.2. 
This  correction  to  the  eigenvalues  are  applied  when  the  absolute  value  of  the  eigenvalues  are 
in  the  neighborhood  of  zero.  The  correction  is  started  from  0.2  as  a  compromise  between  the 
numerical  stability  and  the  accuracy.  The  modification  is  defined  as  a  quadratic  function,  i.e., 

IA'1  =  AIAI2  +  B  (5.35) 
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with  the  constraints 


tt'l  =  0.2 


d  tt'l\ 

^  M  /  A  =  0.2 


=  1 


when  IAI  =  0.2 


This  results  in  A  =  2.5  and  J3  =  0.1. 

5.4  Viscous  Flux  Calculations 

The  viscous  fluxes  are  treated  as  source  terms  and  are  evaluated  explicitly.  The  vis¬ 
cous  flux  terms,  that  involve  the  components  of  the  viscous  stress  tensor,  are  defined  in  Chap¬ 
ter  II.  The  viscous  stress  tensor  in  turn  is  a  function  of  velocity  and  temperature  gradients. 

The  gradients  of  velocity  and  temperature  are  estimated  by  the  application  of  Green’s 
theorem  to  a  control  volume.  In  this  case  the  control  volume  is  taken  as  the  cell  itself.  The 
properties  at  the  nodes  are  computed  using  the  weighted  averaging  procedure  explained  in 
Section  5.3.2. 

The  viscous  stresses  at  the  cell  face  are  taken  as  the  average  of  those  in  the  cells  on 
both  sides  of  the  face.  The  averaging  is  a  good  approximation  for  the  viscous  stresses,  since 
they  are  diffusive  in  nature.  This  gives  a  second  order  approximation  but  the  truncation  error, 
for  the  structured  part,  is  proportional  to  four  time  the  square  of  the  grid  spacing.  The  con¬ 
tribution  of  the  viscous  fluxes  to  the  net  flux  crossing  the  control  surface  are  estimated  as 
shown  for  the  convective  terms  given  in  Section  5.3.1. 

5.5  Time  Discretization 

The  fully  discretized  form  of  equation  (5.10)  can  be  written  as, 

Qn  +  \  _  Qn  k 

Vi  —  —  ^  E-ij  ‘  Hjdsj  (5.36) 

j  =  i 
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The  time  level  at  which  the  right  hand  side  of  equation  (5.36)  is  computed  results 
in  different  numerical  schemes.  The  evaluation  of  the  right  hand  side  of  equation 
(5.36)  at  nth  time  level  results  in  an  explicit  scheme  and  at  (n+l)th  time  level  will 
result  in  an  implicit  scheme.  For  the  explicit  scheme  the  time  step  is  restricted  by 
the  stability  criteria  and  for  the  implicit  scheme  the  flux  vector  has  to  be  linea¬ 
rized  before  the  calculation  of  the  flux  crossing  the  cell  faces.  This  results  in  set 
of  linear  equations  that  have  to  be  solved  to  obtain  the  conserved  variables  at  the 
(n+l)th  time  level.  The  detailed  description  of  these  methods  are  given  in  the  fol¬ 
lowing  subsections. 


5.5.1  Explicit  Scheme 

For  explicit  time  integration,  a  four  stage  Runge-Kutta  method  is  employed  (Jameson 
et  al.  [51]).  The  flow  variables  at  (n+1  )th  time  step  are  obtained  from  the  variables  at  the  nth 
time  step  in  four  stages.  The  time  integration  can  be  written  as 

g(0)  _  q(N) 

e<»  =  e<°>  +  R<e<0)) 

v  i 

e(2)  =  e<  o'  +  a2^R(e(1)) 

(5.37) 

e<3>  =  fi<0)  +  a3^R(e<2)) 

*  i 

q(4)  =  0(0)  +  a4^R(Q(3)) 

*1 

0(^+1)  _  0(4) 

where  R(Q)  is  the  right  hand  side  of  equation  (3.29)  and  the  coefficients  used  are 
aj  =  0.15,  a.2  =  0.3275,  <23  =  0.57  and  04  =  1.0.  These  weighting  coefficients  are 
available  in  the  literature  and  have  been  experimentally  determined  for  struc¬ 
tured  upwind  codes  ( Turkel  and  Van  Leer  [53] ).  For  steady  state  problems  faster 
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convergence  is  achieved  using  local  time  stepping  in  which  the  maximum  permis¬ 
sible  time  step  is  used  for  each  individual  cells. 


5.5.2  Implicit  Scheme 

In  the  case  of  implicit  schemes  the  numerical  flux  crossing  the  cell  face  is  a  function  of 
the  conserved  variables  at  the  (n+1  )th  time  step  (Koomullil  [52]).  For  the  moving  grids  the 
cell  volume  is  changing  at  every  time  step  and  has  to  satisfy  the  geometric  conservation  law. 
The  detailed  description  about  the  problems  associated  with  the  moving  body  problems  are 
given  in  Chapter  V.  The  discretized  equations  then  become 

AOn  k 

Vi-r1-  =  -  Y  n*l.n,ds,  =  R?+1  (5.38) 

1  At  !/  j  j  1 

7=1 

A  linearization  of  the  RHS  about  the  nth  time  level  results  in 


=  R-  + 


'dR? 

dQj 


AQn 


Substituting  equation  (5.39)  into  equation  (5.38),  yields 


(5.39) 


(5.40) 


Using  the  the  summation  convention  over  the  edges  the  above  system  can  be  written 
as, 

t,iaq"  +t(§)  aq" = r‘  (B-4i) 

Utilizing  the  fact  that  $F y  is  a  function  of  conserved  variables  on  left  and  right 
sides  of  each  face,  the  above  system  of  equations  is  written  as, 
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v -  (  dSA 

f',A®  +1  bo:  AQ"<-  + 


wJ AQ*  - * 


(5.42) 


According  to  the  sign  convention  of  the  surface  normal,  the  surface  normal  points 
from  left  cell  to  the  right  cell,  the  cell  number  on  the  left  hand  side  L  and  the  cell 
number  i  are  the  same.  Then  the  above  equation  reduces  to, 


y.  *  /  d$u\  *  /  dTl7\ 


(5.43) 


The  Jacobian  matrices 


can  be  estimated  either  by  an  approximate 


analytic  method  or  by  a  numerical  approach.  These  are  discussed  in  detail  in  Sections  3.4.2.2 


and  3.4.2.3. 


5.6  Newton  Iterations 

For  unsteady  calculations,  it  is  important  to  drive  the  unsteady  residual  to  zero  for  bet¬ 
ter  resolution  of  the  physical  phenomenon.  This  can  be  done  by  using  Newton  iterations.  The 
basic  principle  behind  Newton’s  method  (Whitfield  and  Taylor  [54],  Whitfield  [56])  and  its 
implementation  are  given  below. 

The  implicit  form  of  the  discretized  governing  equations  are  written  as 


%  =  V, 


.p^)+iw+i>=° 


(5.44) 


The  conserved  variable  vector  Qn+ ^  has  to  satisfy  equation  (5.44)  to  achieve  a  time 


accurate  solution.  This  system  can  be  written  as  3&(Qn  )  ~  0.  Newton’s  method  for  this 

system  is  given  by 


36'  (  <2n  +  l’m)(  Qn+hm  +  l  _  Qn  +  \,m  j  _  -36(g"+1’m) 


(5.45) 
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The  prime  denotes  the  derivative  with  respect  to  the  conserved  variable  vector  Q. 
The  above  equation  can  be  expanded  as 


d%  Ann  + 1  ,  m  d3€  Artn+ 1  ,  m  nr  /  nn  + 1,  m  \ 

dQi  Qi  dQj  Qj  -  } 

where  AQn+hm  =  <2”+1,m+1  -  Qn  +  l>m 


(5.46) 

(5.47) 


where  j  represents  the  cells  the  surrounds  the  cell  i.  The  superscript  m  in  equa¬ 
tions  (5.45M5.47)  represents  the  index  for  Newton  iterations.  Using  equation 
(5.44),  equation  (5.46)  is  written  as, 


to 


'Qn+hm  _  r 

)n  i 

k 

Vi 

At 

i 

+  Y  VS"*1'”) 

> 

7=1 

(5.48) 


At  the  start  of  Newton  iterations,  Qn+1, 0  is  set  as  Qn  and  the  linear  system  will  reduce 


r 

t/  k 

£'+XI 

r»#Vl 

,  dQi  / 

*  /ag 
+  2  [h 

„\"  k 

r  A%  =  -  1^(2") 

(5.49) 

7=1 

7  J 

i- 1  V 

'7  7=1 

This  is  the  same  as  the  implicit  scheme  which  is  used  for  the  steady  state  calcula¬ 
tions  as  given  in  equation  (5.43).  As  more  Newton  iterations  are  used,  then 
AQn  +  i,m  approaches  zero  resulting  in  time  accurate  solutions.  For  more  details 
about  Newton  iterations  (Whitfield  and  Taylor  [54]). 


5.7  Approximate  Analytic  Jacobians 
The  flux  crossing  the  cell  face  is  written  as, 
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(5.50) 


5,j  =  ^(Qr)  +  SHGz)  -  I A  I  (Qr  -  Gi)) 


The  Roe  averaged  matrix  is  IAI  is  a  nonlinear  function  of  the  conserved  variables 
at  left  and  right  sides  of  the  cell  face.  For  simplicity  the  approximate  analytic 
Jacobians  are  calculated  assuming  the  Roe  averaged  matrix  IAI  as  constant 
(Whitfiled  [55]).  The  above  equation  can  be  differentiated  analytically  with  re¬ 
spect  to  Qr  and  Qi  resulting 


=  i  /  mhd 

#Qr  2  [  dQR 

dFij  =  i  (vmD 
dQL  2  ^  dQL 


+  IAI 

-  IAI 


2  ( A/j  +  IAI) 
2  (al  ~  IAI) 


(5.51) 

(5.52) 


Substitution  of  the  above  expressions  into  equation  (5.43)  will  result  in  a  block 
sparse  matrix. 


5,8  Numerical  Jacobians 

The  Jacobians  can  be  calculated  from  first  principles.  Each  element  in  the  Jacobian 

&HQ) 


can  wntten  as  a 


v 


dQi 


The  first  order  form  of  these  quantities  can  be  evaluated  using  a 


finite  difference  formula  (Whitfiled  [55],  Vanden  [57]). 

^(<2  +  hep  -  Tj(0 


aij 


h 


(5.53) 


where  ej  is  thej^  unit  vector  and  h  is  taken  as  the  square  root  of  machine  zero. 
The  choice  of  approximate  analytical  or  numerical  Jacobians  has  a  strong  influ¬ 
ence  on  convergence  to  steady  state 
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5.9  Sparse  Matrix  System 


Equation  (5.43)  will  result  in  a  linear  sparse  block  matrix  system  of  the  form 
Ax  =  b.  The  number  of  non-zero  entries  in  each  row  of  A  depends  on  the  number  of  neigh¬ 
bors  of  the  cell,  whose  cell  number  is  the  row  under  consideration.  If  the  cell  i  has  p  neigh¬ 
bors,  then  the  ith  row  of  A  will  have  p+1  non-zero  elements.  The  matrix  structure  with  a 
sample  grid  is  explained  in  detail  in  NO  TAG.  The  symbol  •  in  the  matrix  represents  a  non¬ 
zero  entry.  Each  non-zero  entry  in  this  matrix  will  be  either  a  4x4  block  for  two  dimensional 
cases  or  5x5  for  three  dimensional  cases.  The  resulting  sparse  matrix  is  solved  using  the  Gen¬ 
eralized  Minimal  RESidual  (GMRES)  (  Saad  and  Schultz  [58] )  method.  An  incomplete  LU 
decomposition  is  used  as  the  preconditioner. 


1 
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Figure  5.7  Sparse  Matrix  Structure  Resulting  from  Implicit  Scheme 

5.10  Time  Step  Calculation 

The  time  steps  used  for  the  computations  of  three  dimensional  problems  are  based  on  a 
direct  extension  of  the  stability  analysis  of  two  dimensional  cases.  The  time  step  is  calculated 
using  the  relation  ( Frink,  et  al.  [48] ) 


Matrix  Structure 
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(5.54) 


where 


At‘  -  CFLA ,  +  Bi  +  C, 


A,  =  ( la,l  +  c, )  Sf> 

B,  =  ( lv,.|  +  c, )  Sf> 

C,  =  ( lw,l  +  c, )  S» 


Vi  is  the  cell  volume,  cj  the  local  speed  of  sound,  and  5, and  are  the 
projected  areas  of  the  cell  i  along  x,  y,  and  z  coordinate  directions  respectively. 


5.11  Boundary  Conditions 

The  boundary  conditions  are  derived  from  the  sign  of  the  eigenvalues  and  the  principle 
of  wave  propagation  (Whitfield  and  Janus  [60]).  Local  one  dimensional  flow  is  assumed  for 
these  derivations.  The  boundary  conditions  for  different  flow  situations  are  discussed  below. 

5.11.1  Supersonic  Inflow  or  Outflow 

All  the  eigenvalues  have  the  same  sign  for  the  supersonic  inflow  or  outflow  condi¬ 
tions.  Therefore,  flow  variables  are  specified  for  the  inflow  condition  and  are  extrapolated  for 
the  outflow  conditions. 


5.11.2  Subsonic  Inflow 

In  this  case  four  eigenvalues  are  of  the  same  sign  and  the  fifth  one  is  of  the  opposite 
sign.  According  to  the  sign  convention  the  direction  of  the  unit  normal  vector  to  the  control 
surface  is  always  pointing  outward  from  the  control  surface.  Therefore  the  subsonic  inflow 
will  always  be  in  a  direction  opposite  to  the  unit  normal.  This  will  result  in  the  following 
conditions  at  the  boundary. 
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pb  =  \  0»  +  Pi  -  P(fQ[hx{um  -  Uj)  +  Hy(v oo  -  v*)  +  nz(w oo  -  WL)  ] ) 


Pb  ~  P00 


Pb~  P« 


ub  =  u  oo  —  n 


P°°  ~  Pfc 

*  Poco 


Vfr  =  Voo 


-  Pec  -  Pb 

Uy  Poco 


-  p*-  Pb 
wb  =  -  "z  p0c0 . - 

where  the  subscript  °°  refers  to  the  freestream  conditions,  L  refers  to  the  cell  inside  the 
domain  and  b  is  the  point  on  the  computational  boundary. 


5.11.3  Subsonic  Outflow 

Similar  to  the  subsonic  inflow,  in  this  case  also  four  eigenvalues  are  of  the  same  sign 
and  the  fifth  one  is  of  the  opposite  sign.  Therefore  one  characteristic  variable  is  specified  in 
the  form  of  pressure  and  the  other  four  are  determined  from  the  information  within  the  do¬ 
main.  So  the  boundary  conditions  are  determined  as  follows 


Pb  Pspecified 


Pb  =  Pl  + 


Pb~PL 


'0 


ub 


UL 


+  nx 


PLZlb 

Poco 


vb 


VL  +  hy 


PLZlb 

Poco 


wb 


wL  +  nz 


PLZlb 

Poco 
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5.11.4  Impermeable  Wall 


For  the  impermeable  wall,  there  is  no  flow  across  it  and  so  the  first  three  eigenvalues 
are  zero.  The  fourth  and  fifth  eigenvalues  have  opposite  signs.  These  will  result  in  the  fol¬ 
lowing  relations  for  the  boundary  values. 


Pb  =  Pl+  Poco  ( n*uL  +  nyvL  +  hzwL  +  ht) 


Pb=  Pl  + 


Pb~  Pl 


'0 


ub  =  UL  -  ”x  (  +  nyvL  +  hzwL  +  ht ) 

vb  =  VL  -  hy  (  hxuL  +  hyvL  +  nzwL  +  nt ) 


wb  =  WL  ~  (  *x «L  +  ">-vL  +  +  "f ) 


where  L  refers  to  the  first  cell  from  the  boundary.  For  the  viscous  calculations  the 
velocity  components  on  the  body  are  set  to  the  local  velocity  of  the  body  and  the  pressure  and 
density  are  calculated  in  the  same  manner  as  that  of  inviscid  case. 


5.12  Turbulence  Modeling 

Viscosity  can  be  viewed  as  consisting  of  two  contributions:  laminar  and  turbulent. 
The  laminar  viscosity,  which  is  a  fluid  property,  is  usually  a  function  of  temperature  and  can 
be  estimated  using  Sutherland’s  formula  (  Anderson  et  al.  [41]).  The  turbulent  viscosity,  on 
the  other  hand,  is  a  function  of  the  flow  and  needs  to  be  evaluated  using  some  model.  In  this 
study  the  turbulent  viscosity  is  estimated  using  the  Spalart— Allmaras  one  equation  model 
(Spalart  and  Allmaras  [35])  and  the  Reynolds  stress  is  modelled  using  the  Boussinisq  hypothe¬ 
sis  (  Warsi  [40] ),  i.e.,  it  is  written  in  terms  of  the  turbulent  viscosity  and  the  gradients  of  the 
flow  variables.  Spalart-Allmaras  turbulence  model  is  a  pointwise  model,  which  makes  the 
application  of  this  model  for  unstructured  or  hybrid  grids  easier.  This  model  solves  a  second 
order  partial  differential  equation  for  the  variable  v  and  the  turbulent  kinematic  viscosity  vt  is 
estimated  from  v  by  multiplying  by  a  damping  function^,; . 
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The  different  terms  appearing  in  equation  (5.55)  can  be  divided  into  convective,  pro¬ 
duction,  diffusion,  and  destruction  terms.  So  the  transport  equation  for  the  turbulent  viscosity 
parameter  is  written  as. 

Time  rate  of  change  of  viscosity  parameter  =  -Convection  +  Production  + 

Diffusion  -  Destruction. 
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where  the  convection  is  given  by  V  •  Vv.  Production  of  the  turbulent  viscosity  is 
due  to  vorticity  and  is  given  by  Cbl  S  v.  The  diffusion  term  is  a  function  of  the 
gradient  of  the  viscosity  parameter  and  is  defined  as 

^(V.((v  +  v)Vv))+^(Vvn2 

Neglecting  the  gradient  of  laminar  viscosity,  the  diffusion  terms  can  be  written 
as, 

1 i^(V.((v  +  v)W))-  ^(•'+«0V2v 

The  destruction  term  for  the  turbulent  viscosity  parameter  is  taken  as  inversely 
proportional  to  the  square  of  the  distance  to  the  solid  wall  and  is  given  by, 

f<o  ( v 


dv 

dt 


ReL  \d 


5.12.1  Integral  Formulation  of  the  Turbulence  Model 
Equation  (5.55)  with  the  simplified  diffusion  terms  can  be  written  as, 

(1  +  cbl) 


V  .V  v  +  CM  Sv  +  _ 

—  bl  ReL  a 


(V  .  ((v  +  v)  Vv) ) 


'bl 


(v  +  v) 


R eL 

Integrating  equation  (5.56)  over  a  control  volume  yields 


ReL  \d 


(5.56) 
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+  —  [  V  •  ((v  +  v)  Vv) 


=  -  F  •  Vrrffi  +  CM  SvdQ 

Jo*  L“  Jc 

+  f  + 

L  JQ 

■  i&L(v  +  v')VWfi 

-  fe  U©*» 


(5.57) 


Using  the  divergence  theorem  for  the  surface  integral,  equation  (5.57)  can  be 
transformed  to, 


(1  +C2 


=  -  F  •  Vvd£  +  CblSvdQ 

JQ  Jq  JQ 

(1  +  Co)  f 

+  -R^"f  (V  +  ^  W  ' 

L  J  dQ 

-  ssfsf  <v  +  v2ird!3 


(v  +  v)  Vv  •  n  ds 


W7  I  Hl)da 


(5.58) 


5.12.2  Numerical  Procedure 


As  in  the  case  of  Navier-Stokes  equations,  equation  (5.58)  is  also  solved  using  a  cell 
centered  scheme.  The  values  of  v,  stored  at  the  cell  center,  are  assumed  to  be  cell  averaged 
values.  The  discretized  form  of  equation  (5.58)  is 


v*+1  -  v" 


-Vi  =  -  Y<P+*i  +  U-VMj))dsj+  CbI  (Sv),Vt 
;= i 

(1  +  Ch2)  X 

+  -vwrl(v  +  v)«<-w 
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a  ReL 


(v  +  v)t  V2vt  Vt 


(5.59) 


where  k  is  the  number  of  edges  of  cell  i,  n  is  the  unit  normal  to  the  face,  dsj  is  the 
length  of  the  edge  j,  n ^  is  the  cell  that  shares  the  jth  edge  of  the  cell  i,  and  Vi  is  the 
cell  volume.  The  variables  U+  and  U~  are  defined  as 


U+  =±(U+\U\) 

U~  =±{U  -\U\) 

where  U  is  the  contravariant  velocity.  The  laminar  and  turbulent  viscosity  at  the 
cell  faces  are  taken  as  the  average  values  of  those  on  either  side  of  the  cell  face. 

The  gradient  of  the  turbulent  viscosity  parameter  v  appear  in  equation  (5.59)  in  the 
surface  integral.  So  the  gradient  is  calculated  on  the  cell  face.  The  gradient  at  the  cell  faces 
are  calculated  using  Green’s  theorem  and  a  weighted  averaging  procedure.  The  value  of  v  at 
the  nodes  are  calculated  by  a  weighted  average  of  v  within  the  cells  surrounding  the  node.  The 
weight  is  taken  as  inverse  of  the  distance  between  the  cell  center  and  the  node. 


V  v  = 


<p  vnds 

1  ha 


(5.60) 


The  gradients  at  the  cell  face  are  calculated  based  on  a  control  volume  surround¬ 
ing  that  face.  The  control  volume  is  taken  as  the  volume  enclosed  by  the  nodes 
connecting  the  face  and  the  cell  centers  as  shown  in  Figure  5.8 
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Figure  5.8  Control  Volume  for  Edge  Connecting  Nodes  N1  and  N2 


5.12.3  Explicit  Method 

For  the  explicit  solution  of  the  equation  (5.59),  a  four  stage  Runge-Kutta  method  is 
utilized.  Denoting  the  right  hand  side  of  equation  (5.59)  as  R(v),  the  four  stage  Runge-Kutta 
method  can  be  written  as, 

^(O)  _  ytl 

vW  =  v"  +  ax  AtRtf®) 
v<2>  =  v*1  +  a2At  R(v^) 

v*3)  =  v”  +  a3At  R(i^2)) 

v<4)  =  v"  +  a4At  R(v*3)) 

^n  + 1)  _  ^4) 

where  a1  =  0.0833  a2  =  0.2069  a3  =  0.4265  a4  =  1.0 

5.12.4  Implicit  Method 

For  the  implicit  scheme  the  turbulent  viscosity  parameter  v  appearing  in  the  right  hand 
side  of  Equation  (5.59)  is  taken  at  (n+1  )th  time  level.  Equation  (5.59)  can  be  written  as, 
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(5.61) 


Yi  W" 

Ati 


C(  v"+1 )  +  P(r,+1)  +  Dl(vn+1)  +  D2(vn+1) 


where  C,  P,  Dj,  and  £>2  are  convective,  production,  diffusion,  and  destruction  terms 
respectively.  Linearization  of  the  convective  term  yields, 

C(  v"+ 1 )  =  C(  r )  +  +  -J£-  0V„W)»  (5.62) 

The  Jacobians  appearing  in  the  Equation  (5.62)  can  be  calculated  analytically 
and  are  given  by, 


| £  =  U+ 

BV: 


dC 


dv 


=  U' 


(5.63) 

(5.64) 


n(f) 


Similarly  linearizing  the  production  term,  one  has, 
J7rt+1\  —  T)(  , rifl  \  1  dP 


(5.65) 


P(vn+l)  =  P(vn)+^(Avir 

The  production  is  not  an  explicit  function  of  v,  and  so  the  Jacobian  is  evaluated 
numerically  as  given  below: 

dP  _  P(y  +  h)  -  P(y)  (5.66) 

dv  h 

where  h  is  the  square  root  of  machine  zero.  Linearization  of  the  diffusion  term 


yields, 


D,(v«+1 )  =  £>,(v”)  +  ( Ary  +  <4vn(jy 

KJVi  n(j) 


(5.67) 


where 


&D,  1  +  ch2  A_  cb2 

-r=r-  =  - 5 - -  {Av  •  «)■■  dS:  -  —5 -  A2V:  V: 

dv;  <rRef  lJ  1  aRe,  1  1 


dDx  1  +  Cb2  /A_ 

-r= —  =  — 7T — -  (Av  •  n)::  ds: 
dvn(j)  O  ReL  ->*  J 


(5.68) 


(5.69) 


Linearizing  the  destruction  terms  yields 
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(5.70) 


D2(vn+l)  =  D2(vn)  +  ^(Avl)n 

The  Jacobian  corresponding  to  the  destruction  term  is  also  evaluated  numerically 
as, 

d£>2  =  +  h)  -  D2{v)  (5.71) 

dv  h 

Substituting  Equations  (5.62),  (5.65),  (5.67),  and  (5.70)  into  (5.61)  and  simplifying 

yield. 


=  C(f")  +  P(v”)  +  D^v")  +  D2(v")  (5.72) 

The  system  of  equations  represented  by  Equation  (5.72)  is  a  sparse  matrix  system  and 
is  solved  using  the  GMRES  solver.  The  GMRES  solver  is  also  used  for  solving  the  matrix 
system  resulting  from  equation  (5.72). 

At  the  farfield  and  the  inflow  boundaries,  the  value  of  v  is  specified  such  that  %  =  5 
(Allmaras  [61])  and  at  the  solid  boundaries  the  turbulent  viscosity  is  set  to  zero. 

5.13  Geometric  Conservation  Law 

The  geometric  conservation  law  adds  a  correction  term  to  the  governing  equations  to 
account  for  the  grid  motion  (Thomas  and  Lombard  [62],  Janus  [63]).  This  preserves  the  uni¬ 
form  conditions  when  the  grid  moves,  thereby  avoiding  the  creation  of  spurious  sources  and 
sinks  in  the  flow  field.  The  geometric  conservation  law  can  be  derived  from  the  continuity 
equation  as  described  below. 

Consider  a  control  volume  Q  with  control  surface  dQ,  that  moves  with  speed  c  with 
respect  to  a  stationary  inertial  frame.  Then  the  continuity  equation  can  be  written  in  integral 
form  as  (Warsi  [40]), 
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d_ 

dt 


Q 


i' 

JdQ 


p  dQ  +  (p  j o  (V  —  c) .  nds  =  0 


(5.73) 


Assuming  uniform  conditions  everywhere  in  the  field,  the  above  equation  reduces  to 
dQ  =  J>  c.nds  (5.74) 


d_ 

dt 


JQ 


dQ 


Equation  (5.74)  represents  the  geometrical  conservation  law  in  integral  form. 
This  relates  the  rate  of  change  of  the  control  volume  to  the  orientation  and  veloci¬ 
ties  of  the  cell  faces. 

For  dynamic  grids,  a  term,  accounting  for  the  changes  in  the  control  volume,  appears 
in  the  numerical  discretization  of  governing  equations.  With  this  extra  term,  the  semi-discre- 
tized  form  of  the  governing  equations  can  be  written  as. 


v"A-§ +  osij  *  + 1 


t  I  dfi  +  I  E(Q).nds  =  0 
Q  J  8Q 


(5.75) 


The  second  term  in  Equation  (5.75)  is  calculated  using  the  geometric  conservation  law, 
and  is  evaluated  as 


Q”  '  nds  =  Q"  ]T( xtnx  +  ytny  +  ztnz)dSj  (5.76) 

J dQ  j= 1 

The  above  summation  is  carried  over  the  cell  faces.  Since  the  data  structure  is 
based  on  the  edges,  the  above  integration  is  done  globally,  as  in  the  case  of  the 
flux  summation. 


5.14  Grid  Movement 

For  flow  simulation  over  moving  bodies,  the  grid  has  to  be  regenerated  globally  at 
each  time  step  or  the  grid  points  have  to  be  moved  appropriately  to  retain  a  body  conforming 
grid  where  grid  lines  do  not  cross.  Regeneration  of  grids  for  each  time  step  is  expensive  and 
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interpolation  has  to  be  done  to  transfer  the  data  between  the  grids.  The  interpolation  is  usually 
non-conservative  and  it  reduces  the  accuracy  of  the  flow  simulation.  In  the  later  approach, 
the  grid  points  are  moved  using  different  methods,  for  example  tension  spring  analogy,  poten¬ 
tial  flow  analogy,  etc.  In  the  present  study,  the  tension  spring  analogy  is  used  for  grid  point 
movement  (Singh  et  al.  [64]). 

In  tension  spring  analogy,  each  edge  of  the  grid  is  assumed  to  be  a  tension  spring.  The 
movement  of  the  interior  points  are  calculated  by  solving  the  system  of  tension  springs  when 
the  boundary  points  are  disturbed.  The  motion  of  the  boundary  points  is  computed  either  from 
the  external  forces  acting  on  them  due  to  the  fluid  flow  or  from  a  prescribed  motion.  The 
spring  stiffness  is  taken  as  inversely  proportional  to  square  of  the  length  of  the  edge.  The 
spring  analogy,  by  specifying  the  boundary  displacement,  will  result  in  a  linear  system  for  the 
interior  point  displacements.  This  system  is  written  as 

X  Kij(Axi~Axj)  =  0  (5.77) 

j 

X  Kij(Ayi~Ayj)  =  0  (5.78) 

i 

where  Ky  is  the  spring  stiffness  corresponding  to  the  edge  connecting  nodes  i  and 
j,  and  is  defined  as  ly~2,  where  ly  is  the  length  of  the  edge.  The  summation  varies 
over  all  nodes  that  are  connected  to  the  node  i.  The  sparse  matrix  system  result¬ 
ing  from  the  Equations  (5.77)  and  (5.78)  is  solved  using  GMRES,  a  sparse  matrix 
solver. 

For  the  situations  involving  large  relative  motion  of  the  bodies,  the  quality  of  the  grid 
will  degrade  after  a  few  hundred  time  steps.  At  this  stage  a  new  grid  is  generated  and  the 
solution  is  transferred  to  the  new  grid.  A  first  order  interpolation  is  used  to  transfer  the  solu¬ 
tion  from  old  grid  to  the  new  one.  The  time  marching  is  restarted  with  the  interpolated  solu- 
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tions  on  the  new  grid.  Example  of  the  flow  simulation  over  dynamically  moving  bodies  and 
the  capability  to  handle  the  grid  movement  are  demonstrated  in  Section  6.5. 

5.15  Rigid  Body  Movement 

The  trajectory  of  moving  bodies  is  determined  by  the  laws  of  classical  mechanics.  It  is 
assumed  that  all  moving  bodies  have  six  degrees  of  freedom.  These  six  degrees  of  freedom 
are  the  translation  along  the  three  coordinate  axes  and  the  rotation  about  the  three  coordinate 
axes.  The  movement  of  the  body  is  determined  based  on  the  pressure  distribution  over  the 
body  due  to  the  fluid  flow  and  the  gravitational  force  due  to  the  weight  of  the  body  (Lohner 
[1],  Koomullil  et  al.  [65]). 

The  notations  used  for  rigid  body  motion  are  shown  in  Figure  5.9.  The  position  vector 
of  any  point  on  the  surface  of  the  body  is  expressed  as  the  vector  sum  of  the  position  vector  of 
the  center  of  gravity  and  the  vector  from  center  of  gravity  to  the  corresponding  point 
(Figure  5.9)  and  is  written  as, 

L  =  Lc  +  Lo 

Now  the  velocity  at  the  point  P  is  given  by, 

L  =  Lc  +  i^  =  Vc  +  MX  rQ 

where  r  represents  the  position  vectors  shown  in  Figure  5.9.  and  the  •  represents 
its  derivative  with  respect  to  time.  The  linear  velocity  and  the  angular  velocity  of 
the  body  are  Yc  and  (o_  respectively. 
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Figure  5.9  Notations  Used  for  Rigid  Body  Dynamics 


The  balance  of  forces  and  moments  acting  on  the  body  results  in  the  relations  for  the 
translational  and  rotational  accelerations.  These  expressions  are  written  as 


mVc  =  ]>> 

'  =  mg - 

<p  pnds 

(5.79) 

J  dQ 

0  •  (O  +  f 

Ifa)  X  (o)dn  =  y  r0  X  F  = 

—  <p  p  rQ  X  nds  (5.80) 

Ji 2 

JdQ 

Iyy  +  I & 

1 

1 

r 

where  0  = 

I  xx  "1"  ^zz  lyz 

hj  =  ri  r{p  dQ 

~~  Ixz 

i 

a' 

+ 

JQ 

and  d)  is  the  angular  acceleration  vector. 

For  two  dimensional  cases,  the  second  term  on  the  left  hand  side  of  equation  (5.80) 
vanishes  and  5)  has  only  one  component,  which  in  turn  simplifies  the  equation  greatly.  For 
two  dimensional  cases,  equation  (5.80)  reduces  to 
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As  the  body  moves  under  the  action  of  the  aerodynamic  and  body  forces,  the  second 
moment  of  inertia  changes.  These  changes  are  only  due  to  the  rotation  of  the  body,  as  the 
body  is  assumed  to  be  rigid.  Figure  5.10  shows  the  co-ordinate  axes  before  and  after  the  rota¬ 
tion  of  the  body  through  an  angle  #  with  respect  to  the  center  of  gravity.  The  moments  of 
inertia  with  respect  to  the  new  co-ordinate  axes  are  related  to  that  with  respect  to  the  old  co¬ 
ordinate  axes  by  (Fletcher  [66]), 


Figure  5.10  Orientation  of  Reference 
Axis  and  Rotated  Axis 

Ix  x<  =  I xx  cos 2  6  +  Iyy  sin2  6  +  21^  sin#  cos  # 

7  ,  »  =  I ^  sin2 6  +  Iyy  cos2 6  —  21^  sin#  cos # 

I x'y*  —  (I xx  ~  Iyy )  sin#  cos#  +  Ixy  (cos2 #  -  sin2#) 

The  moments  of  inertia  are  updated  after  each  time  step. 


(5.82) 

(5.83) 

(5.84) 
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6.  Results  and  Discussions 


6.1  Validation  of  Two  Dimensional  Flows 
Two  dimensional  Euler  calculations  are  validated  by  computing  the  flow  over  the 
NACA0012  airfoil.  The  results  are  compared  with  the  experimented  data.  The  grid  used 
for  this  simulation  consists  of  1961  nodes  and  3332  cells  (Figure  6.11).  The  freestream 
Mach  number  is  0.63  and  the  angle  of  attack  is  2.0  degrees. 

The  results  from  the  first  and  second  order  accurate,  inviscid  calculations  are 
compared  with  the  experimental  data  (AGARD)  in  Figure  6.12.  This  plot  shows  a  notice¬ 
able  improvement  in  accuracy  for  the  second  order  calculations  compared  to  the  first  order 
calculations.  The  results  from  the  second  order  computations  agree  well  with  the  bench¬ 
mark  data. 


Figure  6.11  Hybrid  Grid  Around  NACA0012  Airfoil 


The  convergence  history  associated  with  the  implicit  and  explicit  calculations  for 
the  testcase  described  above  is  given  in  Figure  6.13.  The  implicit  scheme  converges  to  the 
steady  state  within  300  time  steps,  while  the  explicit  scheme  takes  25000  iterations  to  re- 
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X/C 

Figure  6.12  Cp  Distribution  Over  NACA0012  Airfoil  at  Freestream 
Mach  Number  0.63  and  Angle  of  Attack  2.0  Degrees 


duce  the  L2  norm  of  the  residual  to  ten  orders  of  magnitude.  The  L2  residual  is  defined  as 


L2  11=  log 


10 


ncells  4 


1 1 >e,/ 

i=l 7=1 


Numerical  Jacobians  are  used  in  the  implicit  calculations  with  a  CFL  of  50 
and  the  explicit  calculations  are  performed  with  a  CFL  of  1.5.  The  explicit 
scheme  takes  12813.76  seconds  for  25000  iterations  on  a  single  processor  of  an 
SGI  R8000,  while  the  implicit  scheme  takes  1443  seconds  for  300  time  steps. 
Thus  the  time  requirement  for  the  simulation  can  be  reduced  by  an  order  of 
magnitude  by  using  an  implicit  scheme. 

The  laminar  flow  calculations  are  validated  using  the  flow  over  a  flat  plate  at  a 
freestream  Mach  number  of  0.5  and  a  Reynolds  number  of  30000.  The  grid  consists  of 
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Iteration  Number 

Figure  6.13  Convergence  History  for  Explicit  and  Implicit  Schemes  for 
NACA0012  Airfoil  at  Freestream  Mach  Number  0.63  and 
Angle  of  Attack  2.0  Degrees 

2911  nodes  and  2800  rectangles.  The  physical  domain  includes  five  plate  lengths  up¬ 
stream  of  the  leading  edge,  four  plate  lengths  downstream  of  the  trailing  edge,  and  five 
units  above  the  plate.  The  first  point  off  the  plate  is  at  a  distance  of  l.OxlO-4.  The  com¬ 
puted  velocity  profile  is  compared  with  the  Blassius  velocity  profile  (Warsi  [40]),  and  is 
presented  in  Figure  6.14.  The  non-dimensional  t]  axis  is  defined  as. 


The  computed  velocity  profiles  are  plotted  in  Figure  6.14  for  axial  location  0.4, 
0.5,  and  0.6  and  these  profiles  are  independent  of  the  axial  location.  There¬ 
fore,  they  satisfy  the  principle  of  self  similar  boundary  layers  (Warsi  [40]). 

The  turbulent  flow  simulation  is  validated  using  the  flow  over  a  flatplate.  The 
dimensions  of  the  flow  domain  are  identical  to  that  used  for  the  laminar  flow  simulation, 
while  using  a  finer  distribution  of  points  near  the  flat  plate.  The  first  point  of  the  plate  is  at 
a  distance  of  l.OxlO-5.  The  freestream  Mach  number  for  this  simulation  is  taken  as  0.5 
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Figure  6.14  Laminar  Velocity  Profile  Comparison 
With  Theoretical  Data 

and  a  Reynolds  number  of  2xl06  was  used.  The  computed  turbulent  velocity  profile  is 
compared  to  theoretical  results  in  Figure  6.15. 

The  next  geometry  used  for  the  implementation  of  the  turbulence  model  is  the 
standard  NACA0012  airfoil.  The  flow  conditions  used  were  a  freestream  Mach  number  of 
0.799,  an  angle  of  attack  of  2.26  degrees  and  a  Reynolds  number  9.0xl06.  The  grid  used 
for  this  simulation  is  a  structured  grid  of  dimension  290x81  (Figure  6.16)  and  is  converted 
into  the  hybrid  grid  data  structure  before  the  simulation.  A  total  of  200  points  were  given 
on  the  surface  of  the  airfoil.  The  Mach  number  contours  are  shown  in  Figure  6.17. 

The  computed  pressure  coefficient  is  compared  with  the  experimental  data 
(Baldwin  and  Barth  [67])  and  the  results  from  Balwin— Barth  computations 
(Baldwin  and  Barth  [67]).  This  is  demonstrated  in  Figure  6.18.  Slight  devi¬ 
ation  of  the  computed  results  as  compared  to  the  experimental  values  is  real¬ 
ized  in  the  simulation. 


Figure  6.15  Turbulent  Velocity  Profile  for  Reynolds  Number 
2.0E+06  Compared  with  Theoretical  Data 


Figure  6.16  Grid  Used  for  Turbulent  Flow  Over  NACA0012  Airfoil 
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Figure  6.17  Mach  Number  Contours  for  NACA0012  Air¬ 
foil  at  Freestream  Mach  Number  0.799  and 
Angle  of  Attack  2.26  Degrees 


The  memory  requirement  for  double  precision  calculations  is  approximately  595 
words  per  cell  for  implicit  scheme  and  102  words  per  cells  for  explicit  schemes.  Not  much 
effort  has  been  put  towards  the  optimization  of  the  memory  requirements. 


6.2  Effect  of  Limiter  Function  on  Convergence 
The  effect  of  different  limiter  functions  on  the  convergence  of  solutions  to  steady 
state  is  discussed  in  the  following  section.  The  geometry  considered  for  this  study  is  a 
scramjet  engine  (Grimes  [68]).  The  grid  used  for  these  calculations  is  shown  in 
Figure  6.19,  and  consists  of  5515  points  and  10404  triangles.  The  freestream  Mach  num¬ 
ber  and  the  angle  of  attack  for  this  testcase  were  3.5  and  0  degree  respectively.  The  pres¬ 
sure  distribution  inside  the  scramjet  engine,  obtained  from  the  simulation,  is  shown  in 
Figure  6.20. 
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Figure  6.18  Pressure  Coefficient  Comparison  over  NACA0012 
Airfoil  at  Freestream  Mach  Number  0.799  and 
Angle  of  Attack  2.26  Degrees 


Figure  6.19  Unstructured  Grid  for  Scramjet  Like  Geometry. 


Convergence  history  for  this  simulation  with  different  limiter  function  options  is 
plotted  in  Figure  6.21.  The  residual  stalls  after  two  orders  of  magnitude  drop  in  the  case 
of  Barth  limiter.  However,  residual  drops  to  the  machine  zero  with  a  higher  threshold  co¬ 
efficient  K  for  Venkatakrishnan’s  limiter.  The  threshold  factor  K  has  a  strong  influence  on 
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convergence,  as  can  be  seen  in  Figure  6.21.  For  larger  values  of  the  coefficient  K,  the 
limiter  function  is  limited  to  places  with  high  gradients.  By  setting  the  threshold  coeffi¬ 
cient  K  to  zero,  the  limiter  is  active  everywhere.  With  a  proper  choice  of  the  threshold 
factor,  Venkatakrishnan’s  limiter  gives  a  better  convergence  than  Barth’s  limiter. 


0.053  0.217  0.380  0.544 


Figure  6.20  Pressure  Distribution  in  Scramjet  Engine  for 
Freestream  Mach  Number  3.5 

The  computed  skin  friction  distribution  over  the  flat  plat  for  a  Reynolds  number  of 
2xl06  and  a  freestream  Mach  number  of  0.5  is  compared  with  the  theoretical  values  in 
Figure  6.22.  The  theoretical  values  are  calculated  using  the  following  expression  (Warsi 
[40]). 

The  Venkatakrishnan’s  limiter  under  predicts  the  the  skin  friction,  but  the 
Barth  limiter  gives  a  non-smooth  skin  friction  coefficients.  For  very  high  val¬ 
ues  of  the  threshold  parameter  (K=50)  the  distribution  of  the  skin  friction 
deviates  from  theoretical  ones. 

The  effect  of  the  limiter  function  on  the  skin  friction  coefficient  is  studied  using 
the  flow  over  the  NACA0012  airfoil,  which  is  used  for  the  validation  of  the  turbulence 
flow  calculation.  The  skin  friction  distribution  over  the  airfoil  is  plotted  in  Figure  6.23. 
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Figure  6.21  Effect  of  Limiter  on  Convergence  to  Steady 

State  Solution 

cy  =  0.02628  Re"1/7 

The  skin  friction  distribution  does  not  change  with  different  threshold  parameters  for  Ven- 
katakrishnan’s  limiter.  But  Barth  limiter  gives  higher  skin  friction  at  the  leading  edge  of 
the  airfoil  and  lower  skin  friction  after  the  leading  edge.  After  the  shock  induced  separa¬ 
tion  at  the  upper  surface  of  the  airfoil  both  the  limiters  gives  the  same  skin  friction  except 
for  K=20. 

6.3  Effect  of  Approximate  and  Numerical  Jacobians  On  Convergence 
The  scramjet  engine  geometry,  shown  in  Figure  6.19,  is  used  to  study  the  effect  of 
approximate  analytical  and  numerical  Jacobians  on  the  convergence  of  solution  to  steady 
state.  The  freestream  Mach  number  of  3.5  and  an  angle  of  attack  of  0.0  degree  were  used 
for  this  simulation.  The  convergence  history  for  this  testcase  with  different  Jacobians  is 
plotted  in  Figure  6.24  and  Figure  6.25.  Figure  6.24  shows  the  convergence  history  for  a 
CFL  number  of  5.  It  can  be  seen  that  there  is  not  much  advantage  in  using  the  numerical 
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Figure  6.22  Skin  Friction  Distribution  Over  a  Flat  Plate  at 
Reynolds  Number  2xl06  and  Freestream 
Mach  Number  0.5 

Jacobian  compared  to  the  approximate  Jacobian.  When  the  CFL  number  is  increased  to 
20,  the  solution  converges  much  faster  for  the  numerical  Jacobian  compared  to  the  approx¬ 
imate  analytical  Jacobian  as  shown  in  Figure  6.25.  The  numerical  Jacobian  calculations 
take  approximately  5%  more  CPU  time  per  iteration  than  the  approximate  analytical  Jaco¬ 
bian.  But  the  total  time  for  the  simulation  is  less  for  the  numerical  Jacobians  because  of 
the  faster  convergence  to  the  steady  state.  The  memory  requirements  are  the  same  for  both 
the  approaches.  During  the  time  integration  the  Jacobians  are  evaluated  at  each  time  level. 

The  same  behavior  is  noticed  for  the  subsonic  flow  over  the  NACA0012  airfoil. 
The  freestream  Mach  number  and  angle  of  attack  for  this  simulation  are  0.63  and  2.0  de¬ 
grees  respectively.  The  convergence  history  for  the  implicit  scheme  with  approximate 
analytical  and  numerical  Jacobians  are  plotted  in  Figure  6.26.  These  computations  were 
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Figure  6.23  Skin  Friction  Coefficient  Distribution  Over  NACA0012 
Airfoil  at  Freestream  Mach  Number  0.799  and  Angle 
of  Attack  2.26  Degrees 


performed  at  a  CFL  of  50.  From  Figure  6.26,  it  is  evident  that  the  numerical  Jacobian 
gives  a  better  convergence  rate  than  the  approximate  analytical  Jacobian. 

For  the  Navier-Stokes  calculations,  the  numerical  Jacobian  does  not  have  a  signifi¬ 
cant  advantage  over  the  approximate  analytical  Jacobian.  But  the  solution  process  em¬ 
ploying  numerical  Jacobians  is  more  stable  compared  to  that  using  the  approximate  analyt¬ 
ical  Jacobians.  This  is  demonstrated  using  the  convergence  history  for  the  laminar  flow 
over  the  flat  plate  in  Figure  6.27.  The  implicit  scheme  with  approximate  analytical  Jaco¬ 
bian  was  not  stable  at  a  CFL  of  20. 

6  4  Examples  of  Mixed  Element  Type  Grids 
Laminar  flow  simulations  were  carried  out  for  a  four  element  airfoil  at  a  free¬ 
stream  Mach  number  of  0.201  and  an  angle  of  attack  of  0.0  degree.  The  Reynolds  number 


Figure  6.24  Convergence  History  for  Scramjet  Geometry  Using 
Approximate  Analytical  and  Numerical 
Jacobians  for  CFL  5 


Figure  6.25  Convergence  History  for  Scramjet  Geometry  Using 
Approximate  Analytical  and  Numurical  Jacobians 
for  CFL  20 
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Figure  6.26  Convergence  History  for  NACA0012  Airfoil  at  Freestream  Mach 
Number  0.63  and  Angle  of  Attack  2.0  Degrees  for  CFL  50 

based  on  the  chord  length  of  the  main  body  of  the  airfoil  was  taken  as  2.0xl05.  The  grid 
used  for  this  simulation  is  given  in  Figure  6.28.  It  consists  of  10445  nodes  and  15603 
cells,  of  which  5274  cells  have  four  sides.  The  pressure  and  Mach  number  distributions 
for  this  case  are  given  in  Figure  6.29. 

The  capability  of  handling  the  cells  with  an  arbitrary  number  of  sides  in  the  flow 
field  is  demonstrated  using  the  grid  shown  in  Figure  6.30.  This  case  involves  cells  with 
four,  five,  and  six  sides.  It  contains  2345  nodes  and  2536  cells.  Some  of  the  cells  have 
hanging  nodes,  where  only  two  edges  are  connected  to  that  node.  The  flow  conditions 
taken  were  a  freestream  Mach  number  of  0.30  and  an  angle  of  attack  0.0  degree.  The 
pressure  contours  for  this  simulation  along  with  the  grid  is  shown  in  Figure  6.31. 

6.5  Unsteady  Flow  Simulation 

The  unsteady  calculations  are  validated  using  inviscid  transonic  flow  over 
NACA0012  airfoil  pitching  about  the  quarter  chord.  A  structured  grid,  shown  in 


Figure  6.27  Convergence  History  for  Flat  Plate  Using  Approximate 
Analytical  and  Numerical  Jacobians  at  Freestream 
Mach  Number  0.5  and  Reynolds  Number  30000 


Figure  6.28  Hybrid  Grid  for  Four  Element  Airfoil 


(b)  Mach  Number  Distribution 


Figure  6.29  Pressure  and  Mach  Number  Distributions  for 
a  Four  Element  Airfoil  at  Freestream  Mach 
Number  0.201  and  Angle  of  Attack  0  degree 
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Figure  6.32,  with  6000  cells  and  6211  nodes  is  used  for  this  simulation.  It  is  converted 
into  the  hybrid  grid  data  structure  format  so  that  it  can  be  used  with  the  solver  developed 
in  this  study.  The  unsteady  calculations  were  started  from  a  converged  steady  state  solu¬ 
tion. 

The  movement  of  the  airfoil  is  prescribed  such  that  the  angle  of  attack  varies  sinu¬ 
soidally  according  to  the  following  relation 

a{t)  =  am  +  a0  sin((ot ) 

where  am  is  the  mean  angle  of  attack  and  is  0.016  degrees  and  the  maximum 
deflection  op  =  2.51  degrees  The  NACA0012  airfoil  is  assumed  to  be  pitching 
at  a  reduced  frequency,  k,  of  0.1628  at  a  freestream  Mach  number  of  0.755. 
The  reduced  frequency  is  defined  as, 


where  a>  is  the  frequency  in  radians  per  second,  c  is  the  chord  length,  and  V*,  is 
the  freestream  velocity.  The  computed  lift  history  is  plotted  in  Figure  6.33. 
The  results  are  compared  with  the  experimental  data  of  Landon[69].  The  lift 
history  attains  the  periodic  nature  after  half  a  cycle  of  oscillation. 

The  pressure  over  the  airfoil  at  two  different  angles  of  attack  are  presented  in 
Figure  6.34.  Figure  6.34(a)  gives  the  pressure  distribution  at  an  angle  of  attack  0.48  de¬ 
grees  when  it  is  moving  up.  Figure  6.34(b)  plots  the  pressure  distribution  at  an  angle  of 
attack  0.98  degrees  when  the  airfoil  is  pitching  down.  Even  though  the  angles  of  attack 
are  almost  the  same,  the  shock  location  differs  due  to  the  unsteadiness  of  the  flow. 

The  residual  history  is  plotted  in  Figure  6.35.  The  residual  achieves  a  periodic 
state  after  50  iterations.  This  simulation  was  done  using  a  CFL  of  1000  and  is  kept 
constant  during  the  simulation  and  it  took  50  minutes  of  CPU  time  on  a  single  R4400  SGI 
processor  for  one  period  of  oscillation.  The  coefficients  of  lift  and  drag  generated  during 
the  pitching  of  the  airfoil  along  with  the  angle  of  incidence  are  plotted  in  Figure  6.36. 
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(a)  Overall  View  of  the  Grid 


(b)  Zoomed  View  at  the  Interface  Between  Layers 


Figure  6.30  Generalized  Grid  Containing  More 
Than  Four  Nodes  per  Cell 


The  above  described  pitching  of  the  airfoil  was  also  simulated  on  a  hybrid  grid 
shown  in  Figure  6.11.  The  computed  lift  history  versus  the  angle  of  incidence  is  plotted  in 
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Figure  6.31  The  Contour  Plot  of  Pressure  Distribution  Around 
Circular  Cylinder  at  Freestream  Mach  Number  0.3 


(a)  Complete  Grid 


(b)  Zoomed  View  of  the  Grid 


Figure  6.32  Grid  Used  for  Pitching  Airfoil 


Figure  6.37.  There  is  a  slight  difference  between  the  computed  solutions  using  the  hybrid 
and  structured  grids.  This  is  attributed  to  the  difference  in  the  grid  resolution  between  the 
two  grids.  Figure  6.33-Figure  6.36  are  plotted  using  the  simulation  done  on  the  grid 
shown  in  Figure  6.32. 


6.6  Example  For  Bodies  in  Relative  Motion 
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Capability  to  do  flow  simulations  involving  bodies  in  relative  motion  is  demon¬ 
strated  using  a  missile  separation  from  a  solid  body.  The  geometry  consists  of  a  flat  solid 
surface  and  a  missile.  The  grid  used  for  this  simulation  consists  of  4824  nodes  and  9216 
cells.  The  freestream  Mach  number  is  taken  as  2.0.  The  unsteady  simulation  is  started 
from  a  steady  state  solution.  The  pressure  and  Mach  number  distributions  for  the  steady 
state  condition  are  shown  in  Figure  6.38.  The  trajectory  of  the  missile  is  calculated  by 
using  the  aerodynamic  forces  acting  on  it  and  its  weight.  During  each  time  iteration,  the 
position  of  the  body  changes  and  the  corresponding  changes  in  the  interior  grid  points  are 
calculated  using  the  spring  analogy  described  in  Chapter  V.  The  grids  are  regenerated 
when  the  cells  get  too  skewed  and  the  solution  is  interpolated  to  the  new  grid.  Regenera¬ 
tion  was  carried  out  three  different  times  during  the  solution  process.  The  grid,  pressure 
distribution,  and  Mach  number  distribution  at  two  different  time  levels  are  shown  in 
Figure  6.39  and  Figure  6.40. 

6.7  Validation  For  Three  Dimension  Flow  Simulation 
Flow  simulations  are  performed  over  an  ONERA  M6  wing  at  a  freestream  Mach 
number  of  0.84  and  an  angle  of  attack  of  3.06  degrees.  The  Euler  calculations  are  carried 
out  on  a  structured  grid  of  size  97x25x17.  Two  different  views  of  the  grid  are  shown  in 
Figure  6.41.  A  structured  grid  is  used  for  this  simulation  to  compare  the  results  with 
those  from  an  existing  code  (CFL3D  (Thomas  et  al.  [22])).  The  flow  simulations  are  ac¬ 
complished  on  the  same  grid  using  both  the  CFL3D  code  and  the  present  program  under 
validation.  Then  the  results  of  computations  are  compared  with  the  experimental  data  in 
Figure  6.42.  It  shows  the  Cp  distribution  over  the  upper  and  lower  surfaces  of  the  wing  at 
45%,  65%  and  90%  span.  The  calculated  values  agree  well  with  the  experimental  results 


and  CFL3D  simulation. 
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Figure  6.35  Residual  History  for  Pitching  NACA0012  Airfoil 


Flow  simulation  was  performed  for  a  complete  F-15  fighter  aircraft  with  faired 
inlet  and  nozzle.  The  freestream  conditions  were  a  Mach  number  of  0.9  and  an  angle  of 
attack  of  4.84  degrees.  The  unstructured  grid  was  generated  at  McDonnell  Douglas  Aero¬ 
space,  Saint  Louis.  A  coarser  distribution  of  grid  points  is  used  for  the  fuselage  and  a  fine 
distribution  over  the  wing  surface  (Figure  6.43).  The  grid  consists  of  44,272  nodes  and 
235,363  tetrahedra  with  6404  nodes  and  12,669  faces  on  the  surface. 

A  four  stage  Runge-Kutta  scheme  is  used  for  this  simulation  because  of  the  un¬ 
availability  of  memory  in  the  local  machine  for  the  implicit  scheme.  Four  orders  of  reduc¬ 
tion  in  the  L2  norm  of  the  residual  is  taken  as  the  criteria  for  attaining  the  steady  state 
solution.  The  convergence  history  is  plotted  in  Figure  6.44.  The  pressure  distribution  on 
the  body  surface  together  with  the  grid  on  the  symmetry  plane  is  shown  in  Figure  6.45  and 
the  Mach  number  contours  on  the  entire  aircraft  surface  are  in  Figure  6.46. 
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Figure  6.37  Lift  History  for  Hybrid  and  Structured  Grids 

The  calculated  pressure  distribution  on  the  surface  is  compared  with  the  available 
experimental  data  (  Michal  and  Halt  [70]  ).  The  pressure  distributions  at  36%  and  59% 
span  of  the  wing  are  plotted  together  with  the  experimental  data  in  Figure  6.47. 

6-ft  Preliminary  Work  On  Parallelization 

The  main  drawback  of  the  flow  solvers  using  unstructured  algorithms  is  the  high 
memory  and  CPU  requirements.  Parallelization  of  the  flow  solver  provides  a  means  of 
overcoming  this  disadvantage.  Preliminary  work  has  been  done  towards  the  paralleliza¬ 
tion  of  the  hybrid  flow  solver  developed  during  this  study.  Initial  results  show  the  poten¬ 
tial  benefits  of  parallelization  of  the  flow  solver. 

The  partitioning  of  the  grid  for  parallel  processing  is  done  using  a  public  domain 
software,  METIS  (Karypis  and  Kumar  [71]),  for  graph  partitioning.  The  graph  corre¬ 
sponding  to  the  grid  is  given  as  an  input  to  the  software  and  is  divided  into  a  number  of 
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(a)  Grid  Used 


(b)  Pressure  Distribution 


(c)  Mach  Number  Distribution 

Figure  6.38  Steady  State  Solution  For  Missile  Geometry 
at  Freestream  Mach  Number  2 


(c)  Mach  Number  Distribution 

Figure  6.39  Grid,  Pressure  Distribution  and  Mach  Number 

Distribution  at  t=177.52 


(b)  Grid  Point  Distribution  Over  and  Near  the  Wing 

Figure  6.41  Grid  Used  for  Flow  Simulation  Over 

ONERA  M6  Wing 
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(c)  Pressure  Distribution  at  90%  Span 

Figure  6.42  Cp  Distribution  Over  ONERA  M6  Wing  at  Mach  Number 

0.84  and  Angle  of  Attack  3.06  Degrees 


Figure  6.43  Unstructured  Surface  Grid  Used  for  Flow  Simulations 

Over  F-15  Aircraft 
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(a)  Pressure  Distribution  at  36%  Span 


X/C 


(b)  Pressure  Distribution  at  59%  Span 

Figure  6.47  Cp  Distribution  Over  the  Wing  Surface  for  F— 15  Aircraft 
At  Mach  Number  0.9  and  Angle  of  Attack  4.84  Degrees 
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domains  (user  specified)  with  minimum  number  of  edge  cutting.  The  details  of  the  graph 
partitioning  algorithms  are  described  in  Karypis  and  Kumar  [71]. 

Different  blocks  of  the  partitioned  grid  are  assigned  to  different  processors.  The 
information  across  the  block  boundaries  are  passed  as  explicit  boundary  conditions.  Mes¬ 
sage  Passing  Interface  (MPI  Gropp.  et  al.  [72])  is  used  to  communicate  between  the  pro¬ 
cessors.  Results  are  presented  for  a  three  element  airfoil  at  a  freestream  Mach  number  of 
0.2  and  an  angle  of  attack  of  16.2  degrees.  An  example  of  the  grid  partitioned  into  four 
domains  is  presented  in  Figure  6.48.  The  complete  grid  contains  29830  cells  and  is  di¬ 
vided  into  four  blocks,  out  of  which  two  blocks  contain  7458  cells  and  the  other  two  con¬ 
tain  7457  cells.  The  Cp  distribution  from  the  sequential  and  parallel  solver,  with  four 
blocks,  is  compared  in  Figure  6.49.  It  can  also  be  seen  from  Figure  6.49  that  the  Cp  dis¬ 
tribution  from  the  sequential  and  parallel  solvers  matches  well. 

The  actual  speed  up  of  the  parallel  solver  is  compared  with  the  ideal  one  in 
Figure  6.50.  The  SGI  processor  (ANDY)  gives  a  better  speed  up  as  compared  to  the  SUN 
(superMSPARC)  processor,  because  the  data  transfer  across  the  processor  is  faster  for  the 
SGI  processor. 


(a)  Grid  Around  the  Three  Element  Airfoil 


(b)  Grid  Partitioned  into  Four  Blocks 

Figure  6.48  Grid  for  a  Three  Element  Airfoil  and  its 
Partition  into  Four  Blocks 


119 


Figure  6.49  Cp  Distribution  for  Sequential  and  Parallel  Solvers  for 
Three  Element  Airfoil  at  Freestream  Mach  Number  0.2 
and  Angle  of  Attack  16.2  degrees 


Figure  6.50  Speed  Up  of  Parallel  Solver  With  the 

Ideal  Speed  Up 


7.  University-Industry  Collaboration 

The  broad  based  objective  of  the  present  collaborative  project  is  to  provide  a  fo¬ 
rum  through  which  university  and  industrial  researchers  can  jointly  pursue  research  and  devel¬ 
opment  pertinent  to  hybrid  techniques  in  computational  fluid  dynamic  with  a  concerted  effort 
in  the  development  of  hybrid/generalized  grid  generation  and  flow  solver  technology.  This  col¬ 
laboration  has  worked  out  very  well  for  both  industry  and  the  university.  The  collaboration  has 
been  proven  to  be  very  beneficial  to  participating  students.  One  Ph.D.,  two  masters,  and  two 
undergraduate  students  participated  in  this  program.  During  the  rise  of  this  project,  the  partici¬ 
pating  industry  management  and  researchers  have  gone  through  various  changes  (e.g.  Boeing 
merger,  leaving  of  Dr.  Rob  Rogers  from  the  Teledyne  Grown  Eng.  Corp.,  changing  interests  and 
priorities  at  McDonnell  Douglas  Sites,  etc.)  These  changes  had  an  impact  on  the  progress  of  this 
project.  However,  all  the  targeted  objectives  of  the  project  were  accomplished  after  additional 
ten  months  of  no  cost  extension  of  this  project. 

Grid  generation  and  flow  simulation  capabilities  have  been  transferred  to  all  in¬ 
dustrial  sites.  In  particular: 

•  HYBGD2D,  HYBFL2D-3D  and  MARCH3D  systems  are  now  fully  operational  at 
the  Teledyne  Brown  Engineering  Corporation.  The  students  and  faculty  have  made 
various  trips  to  Teledyne  Brown  Engineering  to  work  side-by-side  with  TBE  person¬ 
nel  in  validating  these  systems  for  TBE  applications  (Dr.  Sonat  Praharaj  is  the  new 
contact  person  at  TBE). 

•  A  completely  unstructured  grid  around  complete  F15  aircraft  was  obtained  from  Boe¬ 
ing,  St.  Louis  (McDonnell  Aircraft  Co,  contact  persons:  Timothy  Gatzke  and  Ray¬ 
mond  Cosner)  and  the  hybrid  solutions  were  compared  for  validation.  The  modular 
of  grid  generation  technology  and  flow  solvers  have  been  transferred  to  Boeing-St 
Louis  and  Boeing,  Huntington  Beach  (McDonnell  Douglas  Space  Systems  division, 
contact  persons:  Dan  Parish  and  Darrin  Fricker). 

•  The  gird  generation  and  hybrid  flow  simulation  capabilities  were  applied  to  problems 
of  interest  at  Lockheed-Martin  at  Marietta  (Contact  person:  David  Vaughna).  Vari¬ 
ous  helpful  suggestions  in  this  development  were  supplied  by  Lockheed  Co. 

•  The  simulation  capability  of  the  code  has  been  compared  with  the  COBALT  system 
(a  CHSSI  system,  under  development  at  Wright  Patterson  AFB,  contact  persons: 

Don  Kinsey  and  Bill  Strang).  These  results  will  be  reported  as  a  joint  publication  at  a 
later  date. 
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